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ABSTRACT 

A method for calculating the effect of viscosity on low-density flow 
through a nozzle is described, and details of a computer program for 
applying the method are presented. The analysis is based on the slender- 
channel equations, with slip boundary conditions at the walls. The solution 
is started upstream of the throat, using asymptotic results for slow viscous 
flow in a converging cone. An implicit finite-difference scheme is then 
used to calculate the pressure, and profiles of velocity and enthalpy at suc- 
cessive stations along the channel. 

The cases presented show the effects of the nozzle geometry, the 
Reynolds number, and the thermal condition of the nozzle wall. The results 
show that specific impulse is improved by a throat whose longitudinal radius 
of curvature is small, and that exit-area ratios as low as ten can be used 
without serious loss of performance. 

It is shown that, at sufficiently low Reynolds numbers and low exit- 
cone angles, there is no solution of the slender-channel equations in which 
the flow can expand to supersonic conditions. Instead, the boundary layer 
closes, and the solution resembles a viscous subsonic pipe flow. The 
implications of this finding on the upstream influence of the exit-plane con- 
ditions and on the limits of validity of the slender-channel equations are 
discussed. 

Readers desiring to make calculations can proceed directly to 
Appendix A, which is a user's manual for the computer program. This 
program has been deposited at the COSMIC Computer Center, University 
of Georgia, Athens, Georgia, 30601. 
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1. INTRODUCTION 


The rocket engines that are used for satellite attitude control are 
often required to produce a thrust less than one pound force, extending in 
some instances to values as low as 10 Ibf. ^ Such low values of the thrust 
dictate the use of low reservoir pressures and small nozzle scale; the 
result is that the flow takes place at very low Reynolds numbers, where 
viscous effects are felt all across the nozzle cross-section. Thus, it is not 
possible to calculate the performance of these devices by the usual inviscid, 
one -dimensional flow approximation, with a small correction due to a thin 
boundary layer. 

The studies described in this report have led to a method of predict- 
ing nozzle performance in most of the Reynolds number range appropriate 
to microthrust rockets. These predictions show good agreement with 
experiment. The method is a numerical one, and a relatively simple com- 
puter program has been developed for applying the method. The only inputs 
required for a calculation are the reservoir conditions, the nozzle geometry, 
the gas properties (such as the specific -heat ratio), and (if heat transfer is 
allowed) the distribution of wall temperature. The program determines the 
mass flow, and then prints the pressure, thrust coefficient, and profiles of 
velocity, density, and static enthalpy at selected stations along the nozzle. 

A search of the literature was made at the beginning of this pro- 
gram. The major results of that search are presented in Section 2. It was 
decided, on the basis of this search, that the theoretical approach which 
offered the greatest potential was a numerical solution of the slender- 
channel equations with slip boundary conditions. These equations are 
derived in Section 3, and a few details of the finite -difference method are 
given in Section 4. Prior to the present study, no one had demonstrated a 
routine numerical method of solving these equations for arbitrary nozzle 
geometry and reservoir conditions. The key to the successful development 
done in this study was to apply, to internal flows, certain advances that 
have been made in recent years in the calculation of viscous axisymmetric 
wakes. In particular, use has been made of implicit finite differences, and 
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inversion of the large matrices which these differences lead to has been 
done by taking advantage of the tridiagonal nature of these matrices. Per- 
haps most important of all, the existence of a singular point in the slender- 
channel equations has been recognized by analogy with the wake problem, 
and analogous means of passing smoothly through the singularity have been 
developed. 

The approximation inherent in the slender-channel equations is 
essentially the same as that made when one assumes that the flow has a 
thick boundary layer, with an inviscid core. The mathematical approach 
commonly used in this approximation is to solve by an iteration process, in 
which the first iteration assumes a purely inviscid flow. The pressure 
distribution of this inviscid flow is then used to calculate a boundary layer, 
whose displacement thickness is subtracted from the geometric cross- 
sectional area, to provide the basis for the next iteration on the inviscid 
core. This iteration method often fails to converge under low-density con- 
ditions, where an inviscid flow is a very poor first approximation. 

In the present approach, this iteration method is not used; instead, 
the pressure distribution is found by satisfying, at each step along the 
channel, an equation which contains information about the singular point. 

By proceeding in this manner, it has been possible to find solutions at 
lower densities than those of previous solutions, and in fact it has been pos- 
sible to identify the density level at which the boundary layer closes and the 
inviscid core disappears. 

The initial profiles that are used to start the solution at near- 
reservoir conditions are described in Section 5. Results from a variety of 
cases are discussed in Section 6. Comparisons with experimental data are 
given first, as a means of validating the theory. These are followed by a 
set of parametric studies, which show that the maximum thrust is achieved 
at lower and lower exit-area ratios as the Reynolds number is decreased, 
and that for sufficiently low Reynolds number the core of the flow cannot 
expand to supersonic conditions in the diverging part of the nozzle. The 
influence of heat transfer is found to be significant, and thrust improvements 
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found with a heated wall suggest that engines from which heat is rejected 
(as for example when a catalyst bed is used) would have better performance 
if the rejected heat were introduced into the divergent walls of the nozzle. 

The present method has application to the design of low-density wind 
tunnels. The computer program described here has in fact been used for 
this purpose at several laboratories in the United States. 

Under continuing support from NASA, an experimental follow-on 
program is being pursued, in which profiles of density in a low-density 
nozzle will be measured by an electron-beam fluorescence probe. This 
work, which is being done by Dr. Dietmar E. Rothe, will provide a very 
significant check on predictions of the theory described here. In addition, 
it can be expected to yield important data on the structure of the flow field 
at the very low-density levels where no supersonic core is predicted by the 
present model. 
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2. LITERATURE SEARCH 


At the beginning of this program, a thorough literature search was 

made. Approximately two hundred references were studied, some of which 

2 3 

were provided by automated information - retrieval systems. ’ The most 

valuable groups of papers were those which presented analytical methods 

4-47 

of predicting viscous flow, experimental measurements of nozzle per- 

, 48-68 , . . .. 1, 69-72 _ .. 

formance, and general reviews of the problem. Two other 

73-79 24 

groups worthy of mention were those describing thrust transients ' 

80 —86 

and experimental techniques for measuring low thrust levels. 


The thrust range of interest in the present investigation was that 
between approximately 10 ^ and 1. 0 lbf. Except perhaps for the lowest 
decade, most of this range lies in the transitional region between continuum 
and free-molecular flow (see, for example, Ref. 1, p. 1160). Conse- 
quently, in carrying out the literature search, primary attention was directed 
toward analytical methods based on a continuum model. The analytical 

t 

methods were essentially of two types: those based on conventional thin- 

boundary-layer theory, and those based on the slender-channel equations. 


In the first type of approach, it is assumed that the flow has an 
inviscid core that can be treated as an isentropic, one -dimensional channel 
flow. The boundary layer that forms on the walls is treated in a variety of 
ways (e. g. , by integral methods or by local application of similarity solu- 
tions), and its growth rate is usually coupled to the rate of expansion of the 
inviscid core. In developing this approach, various authors have used a 
number of further approximations; important among these are the degree to 
which transverse curvature is accounted for, and the manner of specifying 
the initial state of the boundary layer. Many of these treatments have been 
developed to the point of presenting a working computer program, some of 
which (see, for example, Ref. 33) are very well documented. The major 
shortcomings of these methods, when applied to the microthrust rocket 
problem, is that the boundary layer is usually assumed to have zero thick- 
ness at the throat, that slip boundary conditions are generally not used, and 
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that the rate of convergence of the iterations between the boundary-layer 
and core flows becomes slower as the boundary layer thickens. 

The slender-channel equations, on which the second type of 

43 44 

approach is based, were first proposed by Williams, 1 who simplified 
the Navier-Stokes equations on the basis that the nozzle walls diverge very 
slowly. These equations are formally identical with the boundary-layer 
equations, including the full effect of transverse curvature. They differ 
only in being written in cylindrical coordinates aligned with the nozzle axis, 
rather than in surface coordinates. In this approximation the pressure is 
constant across the channel, varying only with distance along the channel. 
Thus, except for the boundary conditions, the slender-channel equations are 
the same as those often used in studies of the laminar axisymmetric wake. 

Most of the existing solutions of the slender-channel equations are 

43-46 42 

those of the self-similar variety, done by Williams and by Adams. 

These solutions have given important qualitative information about low- 
density nozzle flow fields. However, they are of limited value when applied 
to the direct problem, since the conditions for self -similarity require 
nozzles whose shape changes when the Reynolds number is changed. 

In recent years, many nonsimilar boundary-layer problems have 

been treated successfully by numerical methods (see, for example. Refs. 87 

and 88). Numerical solutions of the slender-channel equations have in fact 

38 39 

been presented by Milligan and by Myers; only a few such results are 

available, however, due primarily to limitations in computational facilities. 

In addition, finite -difference solutions of the laminar axisymmetric wake 

gg -92 

problem have been presented by several groups. 

In the light of these observations, it appeared that the approach 
most likely to produce useful engineering results was to develop a finite- 
difference solution of the slender-channel equations, starting with a given 
nozzle geometry, and given reservoir conditions. It was decided to follow 
this approach, using slip boundary conditions at the wall, since many studies 
have shown (see, for example, Ref. 93) that these conditions, with the 
Navier-Stokes equations, lead to valid results well down into the transition- 
flow regime. 
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3. BASIC EQUATIONS 


The starting point for the derivation of the slender-channel equa- 
tions is the Navier-Stokes equations. In cylindrical coordinates, for steady, 
axisymmetric flow, with zero bulk viscosity, these are: 

Continuity : 



Radial momentum: 
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Energy : 


P u2± 
r 2i 


r 


2r 


- _ ^4 = 

2ir 2\r 



(3-4) 


Williams simplified these equations by assuming that the ratio of 
radial to axial velocity components, and the ratio of axial to radial gradients 
were each of the order of the slenderness ratio of the nozzle, i. e. 

• °(lt) - 
/ ^ % 

This assumption appears to be suitable for flows in which shear stresses 
produced by transverse gradients are much larger than those produced by 
axial gradients. The only region that might appear to be excluded from 
this rule is that near the throat, where the axial acceleration is greatest, 
and where boundary layers are generally thinnest. Even in this region, 
however, the ordering used by Williams is acceptable, as can be seen from 
the following argument: the axial velocity increases from the value CL ^ at 
the throat to about 2 the point where A / accorc ^ n g to inviscid, 

one-dimensional theory. ^ Thus an upper limit for the order of 
can be taken as 
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If in addition it is assumed that ^^/dr — O ( f ^ 

Williams' ordering, namely that 


one is led to 
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9 / Br 

If these approximations, including the assumption that 
are now used in Eqs. (3-1) - (3-4), the leading terms are 


o (■&■), 


M 0** } + 57 + J f=° 


(3-5) 


ilA- 

/ 32 
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- _ ^ 4. _L -£- / u.r 

dz r zr \ l ?r / 


(3-6) 


— fr- = o 


(3-7) 
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-h t> 


it = 


( 3 - 8 ) 


_L A / r A 
r l a 


These equations are formally identical with the boundary-layer equations, 
including the radial dependences that account for transverse curvature. 

The boundary conditions applied at the walls are those corresponding 
to first-order velocity slip and temperature jump. The first of these 
relates the velocity component along the wall to its gradient normal to the 
wall. When this component is expressed as a linear combination of U- and 
T/~ and is then simplified according to the slender -channel ordering, the 
result is 


2. - ot u ' 


UL = 
n/ 



ars -&■ 


( 3 - 9 ) 


The corresponding condition for the temperature becomes 


a. - <*- 


— iY—JLfTT ors 

2 & (Y+i) f] OA ( 


( 3 - 10 ) 


In all of the results presented here, the gas has been assumed to 
follow the perfect-gas law, with constant Prandtl number, constant specific 
heat ratio, and with viscosity proportional to a power of the temperature: 
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^ = — — p -A ; s. = 


CDM51 


(3-11) 


-A -P-\ 

> l H J 

Dimensionless Forms 

l 

The coordinates are made dimensionless by the nozzle throat radius: 


2/ r / 

y = / r M > ^ - /r* 


(3-12) 


and the dependent variables are made dimensionless with respect to reser- 
voir conditions: 


a = “/, 


/IS 1 ; ” 



P= t4 , j> = ^> ,© = % 


(3-13) 


In addition, it is useful to normalize the radial coordinate by its wall value 


at each station 


1- 




(3-14) 


In terms of these variables, the equations of motion become 
Continuity: 


P*- /jLb/\ + r [p£\ + Z!r / £a g = 

0 j fzr » 0 


(3-15) 
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Momentum: 


e \ ** 


u/ 2 S- 
31 , 


Izl iZ 

Alt' d-y- 


+ — l 

3 <[«i z ^ 1 ( 3 ^ 


mergy: 


cr IQ + 2® 

d* 


tL <x ZE 

c^X 


( 3 - 16 ) 


J 

P 

8 *1 (Tw Z 

3^ 1 

zlO 

/Bd 

B<C Z ( 



£_ / ^ 20 

*? ( ^ A. 


( 3 - 17 ) 


In these equations, the density has been eliminated by use of the equation of 
state: 

P = I> O ( 3 - 18 ) 

The transformed radial velocity \^J that appears in these equations is 

W = V - dr, Elk ( 3 - 19 ) 

^ ctx 


Note that \aJ is zero both at the wall and on the axis, and that if U/ is 
zero elsewhere, it implies that the local streamline deflection is equal to 
the wall slope, multiplied by the fractional distance to the wall. 

The only parameter that appears in these equations (in addition to 
, cJ , and Q"^ (x) ) is a Reynolds number based on reservoir 

conditions and the throat radius: 


11 




3 - 


\J Z He 

j^o 


(3-20) 


This parameter is convenient to use, since it involves only conditions that 
are known in advance and are accessible to the designer. 

The boundary conditions in the ( v( ) coordinate system are: 
a \-, = 


(3-21) 


2. - oC u / 2. Y ir 


Ctrs -&j 


y-i z 


(O 


v 2 - 


dCC \ 


\l \ V I 

If the wall temperature is prescribed, the thermal boundary condition is 


4 -/ - = 


(3-22) 


! 2.Y ' IT 

cos 

/n\ 

(30) 

)J (Y-\) 2 

b 

(°\-n) | 



.2.-06, 


If the wall is taken to be adiabatic, this second boundary condition is re- 
placed by (see Eq. (3-31) ): 


— /— -h u z 
3^[ V S. 


On the axis, symmetry requires that 


(3-23) 


y[ — -/ 
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It 


(3-24) 
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Conservation Relations 

The relations that express the global conservation of mass, 
momentum, and energy can be written down either from a direct transforma 
tion of the integral forms of the conservation laws, or by manipulation of 
Eqs. (3-15) - (3-17). The continuity equation, when integrated from zero 
to V ^ , becomes an expression for \/J : 




© J 


>7tV 
^ Y( (X 

~~o 


cty. 


p«-„ 


l£-*l 

O 


(3-25) 


In particular, if the upper limit is set equal to one (where VJ — O ), this 
becomes an expression for the total mass flow: 



tlvi M 

Q 


A = 


YVV 


(3-26) 


The parameter f\ is proportional to the discharge coefficient of the 
nozzle: 



(3-27) 


The constant of proportionality between (\ and C D is: 


Y = 

i. i 

1. 2 

1. 3 

1. 4 

1. 667 

A /c b “ 

0. 06698 

0. 09361 

0. 1133 

0. 1294 

0. 1624 


Detailed derivations are given in Appendix B. 
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V3^ ir* 


A convenient definition of the stream function is 



A UCv[ dA. 

I 0 

f 1 LZ 

o e 


(3-28) 


Integration of the momentum equation leads to an expression for the 
variation of the thrust coefficient: 


ctF 

d-y- 


P (%?) 

d-V- 


4* . 
3 



(3-29) 


where 


F - 




(3-30) 


Equation (3-29) states that the maximum thrust occurs at the point where 
the contribution from the pressure forces acting on the wall is balanced by 
the viscous shear stress at the wall, an observation that has been made by 
many authors (see, for example, Ref. 56). 


If the momentum equation is multiplied by Z UC • integrated, and 
added to the integrated energy equation, the result is a statement that the 
flux of total enthalpy is changed only by heat transfer at the wall: 


d_ 

d-y- 



yiDUO ( O + 





J? 






(3-31) 
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The adiabatic wall boundary condition, cited above as Eq. (3-23), is found 
by setting the right-hand side of this expression equal to zero. 

Streamtube Relation 

In the study of isentropic, one -dimensional channel flows, a very 
important equation is the one that relates the pressure gradient to the rate 

94 

of change of cross-sectional area. 

i (% y ) = 

For the present problem, this expression can be generalized to include the 

effects of viscosity. The generalization to viscous, two-dimensional flows 

95 

was given by Weinbaum and Garvine; the further generalization to the case 
of axisymmetric flows is given below. 

The equation of state is first used to eliminate the density from the 
continuity equation; in the resulting expression, the momentum and energy 
equations are used to replace / <?£ and in favor of 

The result is: 


Vr \ a ) yf? M z tit 



"i- 

The axisymmetric formulas have also been presented in a recent paper by 

96 

Garvine and Weinbaum. 
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The left-hand side of this expression is proportional to the radial gradient 
of the streamline direction, and hence to the axial rate of change of cross- 
sectional area of the streamtu.be. If the streamtube area is taken as 
2~irf Ar , it can be shown that 


± JL 

r dr 



i 

2 tt r Ar 


d- 

d-y 


^ 2-xrr A r ^ 


Thus, the first two terms in Eq. (3-32) are exactly the same as those 

familiar from the case of isentropic, one -dimensional channel flows. The 

additional terms represent the contributions to the area change that come 

from shear, net heat conduction to the streamtube, and heat generation 

within the streamtube due to the conversion of kinetic to thermal energy. 

These terms can be checked against the equivalent expressions derived 

97 

from first principles by Shapiro. 


In terms of the dimensionless variables, the streamtube relation 



l - M Z 

/PM 2 




dp 

eC* 
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(lM) 6 ^ PM Z 




a f 3_(^_ 
© L 9 i V R- 



(3-33) 


If this expression is now integrated across the channel, it becomes a rela- 
tion between the rate of change of the total channel cross-sectional area, 
and the integrated effects of pressure gradient, shear, heat conduction, and 
heat generation. It is useful to solve the resulting expression for the 
pressure gradient: 
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d? = 

cLy. 


li'Pr/ 4- 






l ( 31 


(J 
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d_ (AO 20 
21V ^ 2*1 






dv 


M 


— ^ cr w 


/ - M 
M 2 


n ^ 


(3-34) 


In the limit of infinite Reynolds numbers, this expression reduces to the 
familiar inviscid, one -dimensional channel-flow result. At large finite 
Reynolds numbers, the general character of the solutions for various mass 
flows is qualitatively the same as that of the inviscid limit, namely: at suf- 

ficiently low mass flows, the pressure decreases along the converging part 
of the channel, reaches a minimum at or slightly beyond the throat, and 
then rises. For somewhat higher mass flows, the denominator of Eq. (3-34) 
changes sign upstream of the throat. Thus the pressure gradient becomes 
infinite, indicating that the flow has choked. At some intermediate mass 
flow, it is possible to find a saddle point, where the numerator and denomi- 
nator of Eq. (3-34) vanish simultaneously. This intermediate mass flow is 
the only one that allows an expansion to conditions that are "supersonic", in 
the sense that 


/ - M 2 j 

— Ad.i[CO 

J M 1 

o 

Thus for high Reynolds numbers there is a critical mass flow, which allows 
a supersonic flow. Greater mass flows lead to choking, while lower ones 
yield a subsonic flow in which the pressure reaches a minimum near the 
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throat and then rises. The results given in Section 6 illustrate these types 
of behavior. 

However, the results also show that the saddle point moves farther 
downstream as the Reynolds number is lowered, and that at low enough 
Reynolds numbers it does not occur at all within the range of exit -area 
ratios normally of interest in rocket design. For these cases, there is no 
mass flow that leads to supersonic conditions in the average sense indicated 
above. There is a certain mass flow for which choking occurs at the exit 
plane. For higher mass flows, choking occurs within the nozzle; for lower 
mass flows, the pressure may either reach a minimum and then rise, or it 
may decrease monotonically to the exit plane if the pressure drop due to 
friction is great enough to overcome the effect of nozzle expansion. 

It should be noted that in cases where the no -slip condition is applied 

to the velocity at the wall, the integrals in Eq. (3-34) diverge. In order to 

obtain a determinate result, some account must be taken of the fact that in 

Eq. (3-33) the pressure gradient and shear-stress terms, .as well as the 

heat-conduction and heat-generation terms, cancel in pairs in a region near 

98 

the wall. Thus the simplification employed by Lighthill might be used, in 
which the integration is stopped short of the wall. In the present paper, all 
of the cases studied allow slip at the walls, and the integrals can be carried 
out to vfl = 1 without complications. 


In some cases the rate of rise is sufficient to cause reverse flow near the 
wall. When this happens, the solution must be terminated, since the 
slender-channel equations are parabolic, and cannot allow an upstream 
influence. 
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4. NUMERICAL METHOD 


Difference Equations 

There are many ways of representing the present set of nonlinear 
partial differential equations by an equivalent set of difference equations. 
The choices made here were patterned after those that have been successful 
in earlier boundary-layer work. 

The coordinates are expressed as 

* = +- («-•) A* 5 n = ( L '0 M , L « !, 2 / * • 


For all the calculations reported here, /_ was taken as 101, and 

that constant is built into the computer program described in Appendix A. 
Values of the dependent variables at a grid point are represented by the 


notation: 


a (*, n) = u S 


Crank-Nicholson implicit differences were used, to rewrite the momentum 
and energy equations in the following form 


Af. a L 

U-K 


L + l ~ LZ'k 


l~ r L+l r L-l L + ! L — / 

^ Wy - &K + (Xk + i - CZ <+ , 


i cJ f L-h I L.4-1 L - I 

I y. d-P + ( Ok ) ~ &K + &K + I ~ j£g, 

■2-V d-y £> 


uH 

^(Aa) z 


(4-1) 


_ Lf I _ t_ . C-l _ C+-I 

2 (Au)" 


— L L- I 

2 -f- IZ 


19 




P i ir L Q <+* 0* 


V> lZ L 4- 

x ^ Bt"R 


W/< 

©r 

L-l 

~ Ok + 

04 


44 ^ 

L+l 

L- 

1 L+l 

e< 

- e< 

+ ©<+i 


0 * + "' 


4 ^ 


^ (gT - aTXcC - © 4 ') 


4 ^M)' 


L+ 1 L L - 1 L + l .-. *- _ I- - | 

Ok - 2 0 K + (£>ic +• 0«fi “ -2 0<c+i + 0^+i 


.4- 4 R- v( 


4 ^vl) 2 " 

^r-^-p^r; - ^*+ 7 ) 

4 (^vi) 2 - 


(4-2) 


Note that the wall radius and pressure appearing here are evaluated at a 
point halfway across the step: 
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The velocity boundary condition at the wall is represented by: 
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For an adiabatic wall, the thermal boundary condition is taken as 
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For the case where heat transfer is allowed, this is replaced by 
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Note that in these boundary conditions the pressure and wall radius have 
been determined at the end of the step: 
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Solution of the Difference Equations 


The solution of the finite -differ enc e equations can be determined 
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recursively, following the method outlined by Richtmeyer and Morton 
(see especially Sec. 11. 5). This solution presumes that the profiles at the 
beginning of the step, and a value of / X-Y ’ are known. 
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The momentum and energy equations are rewritten in the form: 
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where the matrix coefficients {given in detail in Appendix C) depend only on 
known conditions at K ■ Richtmeyer and Morton show that the solution 
of these equations can be found if it is assumed that the solution follows the 
recursion 




(4-7) 


Values of the EL and matrices on the axis are found by rewriting the 

symmetry conditions there as: 
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Values of the EL and -p matrices at points off the axis can now be found, 
in the following way: Eq. (4-7), with /_ replaced by /_-/ , is used to 

eliminate and from Eq. (4-6), and the result is com- 

pared with the original form of Eq. (4-7). The result is the recursion 
formula (here the exponent -/ denotes matrix inversion): 
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(4-10) 



The explicit formula for the inverse matrix used here is 
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Once the matrix coefficients are known for L - 1.2, ... f L^txy. -1 , 
the solution for UC^ and (S)^h can be found from Eq. (4-7) by starting 
with their known wall values and working down to L. - 1 . The wall values 
are found by writing the boundary conditions (Eqs. (4-3) - (4-5) ) as 
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where 
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and where, for an adiabatic wall, 
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When heat transfer is allowed, these are replaced by 
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Eq. (4-12) is now equated to Eq. (4-7), evaluated at Lywoc ~ I i the 

solution is 
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where the £p| and p’A/ matrices are defined as 
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Iteration Methods for Finding 


^/s> 


This completes the formulas required to calculate the profiles at the 
end of the step, for a given value of • This parameter must be 

chosen so as to satisfy the three global conservation relations given above, 
as well as the integrated streamtube -area relation. Since the two difference 
equations being used are local expressions for the conservation of momentum 
and energy, it might be expected that the best choice would be to use the 
value of d?Py/ that satisfies either the total mass conservation or the 
integrated streamtube -area relation. The computer program uses both of 
these options; it has been found that solutions determined with either option 
also satisfy the total momentum and energy conservation. Furthermore, 
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except in the near vicinity of the saddle point, solutions found with either 
option also satisfy the other option. Near the saddle point, it is more accu- 
rate to use the streamtube relation; in most of the solutions described here, 
the pressure gradient has been determined in such a way as to enforce mass 
conservation up to the saddle point, and such as to satisfy the streamtube 
relation downstream of the saddle point. The choice of which formula to 
use for is made by selecting the input parameter XIPG . The 

mas s -conservation formula is used (IPG = 4) for X^XIPG, and the 
streamtube relation (IPG = 1 ) for X^> XIPG. Downstream of the saddle 
point, the streamtube relation is used. 

The specific procedure used in enforcing mass conservation is as 
follows: a value of d-P / i s chosen, and the difference equations are 
solved for the profiles of (X and Q at station K+\ • These profiles 

are then used in Eq. (3-26) to calculate the pressure at the end of the step: 


PMEV/ = 


A 


X J' 


(4-18) 


-o © /K + l 

where the integration is done by Simpson's rule. This pressure is com- 
pared with the value that would be found by extrapolating from the beginning 
of the step: 


pxne^p = F k 4 — AX 

d* 


(4-19) 


and the value of J p /dx is adjusted until these two pressures are equal. 
The adjustment of is accomplished by an iteration process: at the 

beginning of the step, an iteration counter called ITER is set equal to 1, a 
first trial value of dF’ / d% > called PG1 , is used, and the difference 
between PNEW and PXTRAP is calculated. A second trial value of 
tf^P/ Xy. ’ ca ^ e< ^ PG2 and equal to 1. 01 PG1 , is then used, with ITER = 2 
and again the difference between the two pressures is calculated. A 
straight-line extrapolation is then used to determine a next guess at 
d? ! > called PG , that should make the pressure difference zero: 
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VLP = PNEW - PXTRAP 
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P<5»-PnJz - P ( Pg,/ - ?e,z) 
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(4-20) 


where PN1 and PN2 are the values of PNEW corresponding to the first 
and second iterations. A third calculation is then made, with 
equal to PG , and ITER '= 3 . The results of this third iteration are 
checked to determine whether the absolute value of PNEW - PXTRAP is 

_3 

less than 10 of the pressure at the beginning of the step. If it is not, the 
iteration counter is increased by 1, and a new trial value PG is calculated 
using as base points for the straight line the iteration just completed, and 
the prior iteration that came closest to zeroing the quantity DLP . A max 
imum of ten iterations is allowed; convergence is usually attained on the 
third or fourth iteration. 

The same procedure is followed when is chosen so as to 

satisfy the streamtube -area relation. On the first two iterations, ^/dy 
is set equal to PG1 and PG2 , the right-hand side of Eq. (3-34) is calcu- 
lated, and called PGNP1 : 
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In these formulas, Simpson's rule was used, and derivatives were calcu 
lated by the formulas 
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The integrands for Simpson's rule were calculated as shown in the above 
formulas, except for the shear and conduction integral at the wall, where 
special formulas derived from the differential equations of motion were 
used, namely: 
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It is clear that the averaging done in the above formulas could have been 
carried out in many other ways. No attempt has been made to determine 
the optimum choices. 

After PGNP1 has been determined, the difference between PGNP1 
and the value of used in calculating it is found, and the results of 

the first two calculations are used as the basis for a straight-line estimate 
for the next iteration: 
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The solution for this third iteration is checked to determine whether DPG 
is less than 1(J times the value of at the previous step. If it is 

not, the iteration counter is increased by one, the formula above is used to 
calculate a new trial value of PG , using as base points the trial just com- 
pleted and the prior iteration with the least absolute value of DPG . A 
maximum of ten iterations is allowed. 

In general, the values of that zero the quantities DLP and 

DPG are different. When using either of the formulas, the mismatch in the 
other quantity (DLP/p^ or DPG/ ( e ^^) ) is typically the order of . 

one percent, for values of Ay = 0. 1 . 

Radial Velocity Calculation 

During each iteration, with either of the pressure-gradient 
algorithms, the radial velocity is calculated by the finite -difference approxi- 
mation of Eq. (3-25): 
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The actual values of vj used in advancing the solution from one step to the 
next were taken as the average of these values and those of the previous 
station: 
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(4-33) 
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Until the iteration on d-P / Jby. is completed, these values are temporarily- 

stored in the E12 array. This averaging was done in an effort to be con- 
sistent with the Crank-Nicholson averaging used in finding (X ; an earlier 
version of the computer program which did not use this averaging for Ia/ 
sometimes developed divergent oscillations in for mass flows 

near the critical value. 

At the beginning of the iterations at each step, the initial guess 
(PG1) used for JLP/^y was taken as the previous value of ^ , 

multiplied by the ratio of the two previous values: 
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The pres sure -gradient ratio PGRT is set equal to 1 at the initial station. 


Determination of the Mass Flow 

Once conditions at the initial station are provided, the procedure 
described above can be used to find the solution for a given mass flow, 
given Reynolds number, given nozzle geometry and wall -temperature dis- 
tribution, and given gas properties. As the solution proceeds toward the 
throat, dP becomes more and more negative; usually it reaches a 
minimum value very near to the throat, and subsequently increases (i. e. , 
becomes less negative) downstream of the throat. Following this, either of 
two types of behavior is usually observed: in the first, dP continues 

to increase, and eventually passes through zero, denoting a minimum in the 
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pressure. This corresponds to a change of sign of the numerator of Eq. 
(3-34), and indicates that the mass flow being used is less than critical - 
i. e. , it is associated with a generally subsonic flow. 

In the second type of solution, which occurs at higher mass flows, 
begins to become more and negative, until a station is reached at 
which there is no solution that will match PNEW and PXTRAP or the left 
and right sides of Eq. (3-34). This condition corresponds to a change of 
sign of the denominator of Eq. (3-34), and indicates that the flow has 
choked - i. e. , the negative pressure gradient required for the specified 
mass flow is essentially infinite. 

In general, it would be expected that some intermediate mass flow 
could be found, for which the numerator and denominator change sign at the 
same point. Of course, it is extremely unlikely that the numerical solution 
would pass smoothly through this saddle point, and it is necessary to force 
the solution through, once the mass flow has been found to a suitable accu- 
racy. The means of carrying the solution through the saddle point is 
described later; for the moment, it remains to describe the logic used in 
determining the critical mass flow. 

For this purpose, an iteration on the value of /? is performed; an 
upper and lower bound for the mass flow (called (\ I and AO respectively) 
are specified as input data, and the first solution is calculated for a value 
of f\ halfway between these bounds. When this solution displays one of 
the two types of behavior described above - i. e. , when its mass flow has 
been identified as either subcritical or supercritical, then its value of /) 
is used as a new lower or upper bound, and a new solution is calculated for 
a value of f\ halfway between the (updated) lower and upper bound. This 
iteration process is continued until the difference between the lower and 
upper bounds is less than an input value called ATEST . When this condi- 
tion is met, the indicator IGOTA is set equal to 1, and a final solution 
pass is calculated. 

The symptoms that are used to detect whether a given mass flow is 
subcritical or supercritical are as follows: The mass flow is declared to 


'V 
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1 . 

2 . 

3. 


be subcritical if becomes positive, and supercritical if any of 

these four conditions occur 

more than ten iterations are used in finding 

the integral J ^ falls below 0. 5 

if , having passed its first minimum near the 

throat, turns more negative for 3 consecutive steps 

4. if p falls below 0. 1 upstream of the minimum in 

The third and fourth of these conditions involve the indicator IXSTAR 
which is equal to 1 upstream of the minimum in > and is set equal 

to 2 downstream of this minimum. The test to determine whether IXSTAR 
should be changed from 1 to 2 is not applied until after 5 steps have been 
taken, since occasionally some slight oscillations in d-P/dy. occur near 
the initial station. 

In cases where ten iterations on d-P/fa are exceeded, it usually 
happens that the curve of DLP versus W has no zero: 



In such cases, the trial values assigned by Eq. (4-20) for some of the iter- 
ations may become large positive numbers, as indicated in the sketch 
above. But large positive values of can produce reverse flow near 

the wall, which cannot be allowed, by the parabolic nature of the slender- 
channel equations. To avoid this possibility, the pressure gradient is 

restricted to values less than +1. If any of the trial values exceed this 

-3 

limit, they are reduced by a factor 10 , and are simply allowed to count 

as another iteration. 
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Passage through the Saddle Point 

The procedures described above are used repetitively, until the 
mass -flow parameter f\ has been determined to within the preassigned 
accuracy ATEST . The indicator IGOTA is then set equal to 1, and a 
final pass is made. This integration is interrupted slightly upstream of the 
saddle point, and the pressure -g radient distribution is then chosen so as 
to force the solution to lie along a straight line through the saddle point, in 
the plane of the numerator and denominator of Eq. (3-34). As soon as this 
fitted portion of the solution passes through the saddle point, the iterative 
solution for d.P J/ (satisfying the streamtube relation) is resumed. 

The point at which the solution is interrupted is determined in the 
following way: if the solutions for various mass flows are plotted in the 

plane of DPNUM vs AMM2 (except for a factor , this is essentially 

the plane of the numerator and denominator of the streamtube relation, and 
is referred to hereafter as "the saddle-point plane"), their behavior is 


jpP/s/O'AY 


fi MMZ 



The solution for the final value of f\ will be one or the other of these two 
types. If it is supercritical, it should be interrupted at or near the point 
where supercritical behavior was detected on the previous pass that set the 
upper bound (\ ! . This is accomplished by the quantity XCUT . When- 

ever supercritical behavior is found, XCUT is set equal to the value of X 
one step upstream. On the final solution, the straight-line passage through 
the saddle point is initiated as soon as X. passes XCUT . 
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If the final solution is of the subcritical type, it should be inter- 
rupted when the quantity AMM2 reaches its minimum. It has been observed 
on runs with wall cooling (in particular, case 12 of Section 6) that local 
minima can occur far upstream of the saddle point. To discriminate 
against this possibility, the program as presently written will interpret a 
minimum in AMM2 as indicative of nearness to the saddle point only if it 
occurs for y > - OS • 

When either of these conditions is met, the indicator IPG is set equal 
to 2, and a straight line is fitted between the point just calculated (AMM21 , 
DPNM1 ) and the saddle point (0. 5, 0) 



&PMUM - J>Pr^JM\ - /IMM 2 / 

— (4-35) 

-SPaJ/vu OS - /IMM Zl 

At a given value of X » it is true to a good approximation that DPNUM 
and AMM2 are linear functions of the trial values of < and hence 

of each other. If point A locates the coordinates of the first iteration 
(DPDX = PG1) and point 3 the second (DPDX = PG2), then a straight line 
connecting these points is 

IppfJt'M -DP*Jo>MA MM2 -^MM2^ 

= (4-36) 

DP 6 -~D?rJU MA /)MMZ6 - 
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The intersection of this line with the straight line through the saddle point 
occurs at point C , where the value of AMM2 is 

2 C — 


Z>ptJM l - DPn/l/M/4 +- 
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To the extent that these quantities are linear with dP / jl^. • the value of 

dP/jly. associated with point C would be 


p<$\ • (/iMtfZQ - /5/VlM2c)+p^2.- (AM A A2c - 
DPDY C ~ 


(4-38) 


This value is used for the third iteration, and it has been found that the 
resulting solution lies sufficiently close to the straight line through the 
saddle point that no further iterations are necessary. 

This three -step placement of the solution is continued until it passes 
through the saddle point, after which an iterative solution of the streamtube 
relation is resumed. This could be accomplished by setting IPG equal to 3; 
however, it has been found that the rate of convergence of the iterations at 
this first station downstream of the saddle point is improved if the following 
procedure is followed: the indicator ISW is set equal to 2, and the first two 

iterations (points /} and B ) are calculated. Then IPG is set equal to 3, 
and the value of d-P / c j(^ to be used on the third iteration is found by solv- 
ing the quadratic equation that follows from the (roughly) linear dependence 
of AMM2 and DPNUM on d-P/j^ - i. e. , since 
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where 
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it follows that the streamtube relation can be approximated by 
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The solution of this equation is 
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where 
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This quadratic solution could in fact be used throughout the program, 
in place of the straight-line extrapolation presently incorporated. It is not 
clear how much of a saving in time could be achieved by doing so; however, 
at the first station downstream of the saddle point, there are instances 
where ten iterations would not be enough for convergence of the solution if 
the more nearly correct location of the solution were not determined by the 
quadratic formula very early in the iterations. 
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Once the solution has converged downstream of the saddle point, the 
integration continues until the maximum X -location (XMAX) is passed. 





5. INITIAL CONDITIONS 


In order to start the calculation method described in the previous 
section, it is necessary to specify a value of P , and profiles of iX and 
(C) , at some initial station. In applications of thin-boundary-layer theory 

to the present problem, this initial station is usually taken at the geometric 
throat of the nozzle; there the boundary layer is assumed to have zero thick- 
ness, and the pressure, velocity, and enthalpy are assumed equal to the 
values they would have in an isentropic, inviscid, one -dimensional channel 
flow. 

In the present study, where results are sought at very low Reynolds 
numbers, this approximation is not acceptable; instead, the solution must 
be started far upstream of the throat. Asymptotic solutions of the slender- 
channel equations for slow flow in a converging cone were derived for this 
purpose. Analogous solutions of the Navier -Stokes equations for incom- 
pressible flow have been presented by Ackerberg. 

The initial solution was derived by expanding all of the dependent 
variables in inverse powers of the nozzle radius: 


p= ' + £ - z:k - 
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tx = z. EiM , */= 

^=2 X? hj-z 

The expansion for (X begins with ^ = 2 , since any contributions from 
lower values of A / would not satisfy mass conservation (see Eq. (3-26)). 

If this term is then substituted into the continuity and momentum equations, 
it is found that the first nonzero coefficients in the expansions of P and 
\tJ are 7T^ and Wj ('A.') • The fact that (a - Q implies that the 

flow resembles a conical sink flow in the leading approximation: 
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Finally, substitution of these expressions into the energy equation reveals 
that is the first nonzero coefficient in the enthalpy expansion. 

When the boundary conditions are expanded in powers of • 

the first few velocity terms are found to satisfy: 


(X 2 d) = o , 
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(5-2) 


The first term in the expansion of the enthalpy profile must satisfy either 
T 4 .( 1 ) = o if heat transfer is allowed, or 7^/ ( 7 ) = o if there is no 
heat transfer. 

The ordinary differential equations for LE*/ (~*l) and T^/ (\) can 
be integrated explicitly.’ 1 ' The resulting expressions contain as parameters 
the coefficients ; these latter quantities are found in terms of the 

given mass flow from the expansion of Eq. (3-26), i. e. : 



These equations are presented in 


detail in Appendix D. 


(5-3) 
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The results for the first few terms are 
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Further terms are given in Appendix D. The above result for 
is valid for either thermal condition of the wall. 


Simplified Initial Conditions 

It is desirable to specify the initial conditions to a level of accuracy 
that is commensurate with the accuracy of the finite -difference scheme. 
From an inspection of Eq. (5-3), the indication is that mass conservation 
requires the LC -profile to be specified through order and the 

pressure through order ^ • When this was done in an early version 

of the computer program, however, it was found that the higher-order con- 
tributions were masked by the truncation error in the finite -difference 
results - i. e. , if cLP > J / dy. is specified through order . a 

Simpson 's -rule integration of the resulting UO -profile at x+HX is in 
general different from the right side of Eq. (5-3) by an amount of order 


-r 


greater than ■ For this reason, the iterative method of finding 

^^/dy. ' described i* 1 Section 4, was chosen. This has the effect of doing 

numerically to what the series coefficients TTy would do analyti 

cally. In this method, the series solution is used only to give the starting 
profiles, and an initial guess at ^P/dy • 

To be consistent with this procedure, the starting conditions were 
simplified to 
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This approximate velocity profile fails to conserve mass at the initial sta- 
tion; the error introduced is typically less than one percent. 

For all of the solutions presented in this paper, the initial station 
was chosen as Y e - -5. 
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6. RESULTS 


In the subsections that follow, the nozzle geometry used is 
described first; following this, the general features of the flow fields found 
by changing various parameters are described. Results that relate speci- 
fically to performance characteristics (mass flow, thrust, specific impulse) 
are deferred to the end of this section. 

Geometry 

The particular geometrical configuration used in the present study 
consisted of a convergent cone and a divergent cone, connected by a constant 
radius -of -curvature -section: 



The coordinates are given by 

! 

j X 4 sm : 

! xj r ~ X # + I + £?. ( i - cos - sia & 

\A J l ' \ * / 


w i 


45 


■ V 


\ 

■ v . 


fi 


5ltA 9[ < X < ^ Sin ; 




/ 



(6-1) 


X ^ /?, 1% : 

CT^/ - >C -/i^^ + 1 -hi?, (/- cos 9^ - Sir) % -hast -Q j 

= -Wl 

where 



This geometry is built into the program as a subroutine. Other geometrical 
configurations could be studied without difficulty. 

Cases Studied 

The computer program described earlier has been used to investi- 
gate the effects of varying some of the parameters of the problem. Table 1 
lists the values used for the various parameters. For all cases, the values 
= 1.4, = 0. 75 , u) - 0. 9 , 0(u, - OC r = 1.0 were used. In general, 

a step size of A X = 0. 1 was used; however, for runs with the sharp 
throat ( = 0. 5 ) , it was reduced to A X = 0. 01 in the region 

- 0. 3 ^ X ^ +0. 3. The calculations were performed on an IBM 360, 

Model 65, in double precision. The program listed below is written in 
single precision, and check runs have shown that the results are essentially 
unchanged. A typical case, where the mass flow is unknown to plus or minus 
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20 percent (i. e. , A f\ originally is the order of 0. 04) takes between 5 and 
10 minutes, depending on the step size and the amount of output desired. 

Comparison with Experiment 

Case 1 provides a direct comparison with the experimental data 
6 B 

reported by Yevseyev, who carried out a very detailed probing of the low- 
density flow in a conical nozzle. Results for the centerline Mach number 
(taken from Yevseyev's Fig. 8) are presented in Fig. 1; the present calcu- 
lation is 6 percent high at X = 6 , and about 3 percent low at X - 15. A 
more stringent test of the calculational method is shown in Fig. 2 where the 
measured and calculated velocity profiles at / = 11. 1 (see Yevseyev's 
Fig. 7) show excellent agreement. 

It should be emphasized that the calculations shown contain no 
approximations beyond those implicit in the use of the slender -channel 
equations; they constitute the solution of the direct problem, for a specified 
nozzle geometry and a specified set of reservoir conditions. No attempt has 
been made to improve the agreement by varying the gas -property parameters. 

Effects of Upstream Geometry 

Cases 2 and 3 are identical except for the value of the inlet angle; 
the results of these two runs are indistinguishable for X o , and were 
taken as evidence that the upstream geometry has little effect on the flow. 
Obviously, much more evidence is needed before the quantitative limitations 
of this statement are clear, but it was felt that exploration of this effect 
could be deferred until after the variation of other parameters had been 
studied; all of the other cases described here used = -20° . 

Effect of Reynolds Number 

Cases 2, 4, and 5 show the effect of the Reynolds number. Before 
discussing this effect specifically, it is well to examine how the iteration 
process led to determination of the eigensolution. For case 2, Fig. 3 shows 
the variation of P with X at various values of the mass flow; these 
exhibit the two types of behavior described earlier. This behavior is shown 
more clearly in Fig. 4, which is the saddle -point plane. After the iterations 
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shown had been completed, the value /} = 0. 121 was selected as the eigen- 
value of /■) , and this solution was interrupted at X = 0 ; it was then 

placed on a straight line through the saddle point, as shown in Fig. 4. The 
result of this fitting is shown in the P ,y. plane in Fig. 3. Note that the 
saddle point is crossed about one radius downstream of the throat; there- 
after, the iterative method of finding is resumed, and continues to 

run stably. 

An isometric view of the velocity and enthalpy profiles at various 
positions along the nozzle is shown in Fig. 5. 

Lowering the Reynolds number S> to 800 has the effect of moving 
the saddle point slightly farther downstream (case 4), and a further reduc- 
tion to 3 - 400 moves it out beyond twenty radii (case 5). For this case, 
the solution in the plane of Fig. 4 moves along the straight-line segment 
toward the saddle point, but never reaches it. 

Effect of Exit-Cone Angle 

Cases 4, 6, 7, 8, and 9 show the effect of changing the angle , 

at a constant Reynolds number 3> ~ 800. Reducing this angle inhibits the 
expansion to supersonic conditions, and so has the same effect as lowering 
$ . As becomes smaller, the saddle point moves farther down the 

nozzle, passing beyond twenty radii for ^ 18° . There is obviously a 

range of expansion angles and Reynolds numbers for which there is no unique 
supersonic solution of the problem determined by passage through the 
singular point. 

Extended Study of the Nonsingular -Solution Region 

The conditions of case 10 lie well within the region where there is no 
saddle point. A large number of computer runs was made for this case, and 
the- resulting pressure distributions along the nozzle are shown in Fig. 6. 

At the lowest value of the mass flow ( /) = 0. 04 ) , d.3 / is negative out 
to forty radii; at somewhat higher mass flows, the pressure reaches a 
minimum, and then rises slightly as the flow continues to decelerate. 
Ultimately a pressure maximum is reached, and the final state of the flow 
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is a Poiseuille velocity profile, with the static enthalpy uniformly equal to 
its wall value. 

For sufficiently high mass flows, a choking condition is encountered; 
however, in contrast to the higher Reynolds number, higher cases, 

the locus of zero pressure gradients does not intersect with the locus of 
infinite pressure gradients. The solution for f\ = 0. 08764252 extends all 
the way to forty radii with negative. Figure 6 shows this solu- 

tion passing between the loci of the zero- and infinite -slope points. 

For values of /? greater than about 0. 086, these solutions exhibit 
a zone of supersonic flow near the axis, and slightly downstream of the 
throat. Contours of constant Mach number are shown in Fig. 7 for the case 
f\ - 0. 08764252. A smooth transition from supersonic back to subsonic 
conditions is predicted. This type of transition is possible within the frame- 
work of the slender -channel equations, where the axial gradients in the 
stress tensor (which could give rise to a shock -wave type of transition) have 
been neglected in comparison with transverse shear stresses. The Mach 
number distribution along these streamlines is analogous to that of a 
streamline entering the boundary layer on a flat plate in a supersonic flow; 
the Mach number decays below one as the streamline is carried deeper into 
the slower part of the boundary layer. 

Figure 7 also shows the locus where the axial velocity reaches 99 
percent of its centerline value. This locus shows how the boundary layer 
essentially closes on this flow. In all of these flows, the decreasing pres- 
sure generates two opposing effects: the first, appearing in the inviscid 

terms of the streamtube relation, is to accelerate the flow to a supersonic 
state. The opposing effect, appearing in the viscous terms, is to thicken 
the boundary layer. The conditions of case 10 lie in a region where the 
latter effect dominates the former; no acceleration to a permanent super- 
sonic state is possible within the framework of the slender -channnel 
equations, and the flow undergoes a transition back to subsonic conditions. 

Whether a smooth transition of this type is realized experimentally 
is, to the author's knowledge, an open question. Measurements aimed at 
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resolving this problem are needed; in addition, further analysis in which the 
longitudinal stress contributions are retained would be very useful. 

The precise boundary above which a unique supersonic solution can 
be found has not yet been adequately determined; the present results sug- 
gest that this will be the state of the flow if S> is greater than 800 and -&■ 

2* 

is greater than about 20°. This value of & corresponds to a throat 

Reynolds number around 300; using Sutherland and Maes' rule of thumb (see 

-4 

Ref. 1, p. 1160), this corresponds to a thrust level in the vicinity of 10 lbf 
for a typical nozzle size. 

It is interesting to note that the range of £> below which no super - 

20 

sonic solution was found corresponds roughly with Smetana's estimate 
of the Reynolds number at which the boundary layer fills the throat. 

Assuming that the displacement thickness is 0. 27 times the velocity boundary- 
layer thickness when this occurs, Smetana derives the following formula 
for the Reynolds number at which the boundary layer closes 



U-*. • 2. 

h° 



0.(,s~ 


0.3? 


For the constant , he recommends a value on the order of 4. In 

the present notation, his result is 
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Approximating the density and velocity ratios by their isentropic -flow 
values for j = 1.4 , this becomes 
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Smetana's value of ^ 57*4 = 4 corresponds to 5* CC- 110 , while a value 

3 ^ ~ 9 would recover the present result 3* ^00 • ^ larger value 

of ^ is probably more correct for the present configuration, which has 
a longer entrance than the one used in Smetana's analysis. 

Even if these predictions of an imbedded supersonic zone are accepted 
as having physical significance, there still remains a problem of predicting 
nozzle performance when the pressure imposed at the exit plane is extremely 
low. Figure 8 shows the slender-channel prediction of mass flow as a 
function of exit-plane pressure, if the nozzle is cut off at X = 40. The 
point to be noted is that no solutions were found for which 0 was less 

than around 0. 09. It seems probable that for values of /} between 
0. 08764252 and 0. 08764343, solutions could be found in which choking occurs 
arbitrarily far downstream. However, it may be that some mechanisms 
which are not included in the slender-channel equations must be invoked in 
order to predict how the flow adjusts to the exit conditions encountered in 
the space environment. 


Effect of Exit-Plane Pressure 

The fact that large subsonic regions are predicted suggests that a 

strong upstream influence, of the exit-plane conditions may exist under 

101 

certain conditions. In an earlier publication of the present results, the 

speculation was advanced that this mechanism might explain the anomalous 

j . , 57, 58, 102 

effect of the exit-plane pressure reported by several authors. 

Further evidence has since suggested that this mechanism does not explain 

the observations when the flow has a supersonic core; whether it plays a 

role at the lower Reynolds numbers, where the slender-channel equations 

have no solutions with a supersonic core, remains to be determined. 


The principal item of new evidence is the observation, by the group 

at AVCO, * that experimental arrangements in the test cell can affect the 

-4 

measured thrust if the cell pressure is greater than 10 torr. Values in 


Private communication from R. R. John and Walter Davis, July 8 , 1969. 
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this range have previously been recommended as the maximum desirable 

cell pressure in all three of the references cited above; the more recent 

AVCO observations appear to be the first in which it is clearly shown that 

different arrangements of the same equipment will indicate different thrust 

-4 

levels at cell pressures in excess of 10 torr. 

A second item of new evidence is that calculations made with the 
present computer program for one of the experiments reported in Ref. 102 
show good agreement with the measured performance at low cell pressure. 
The calculations were done for the parameter values listed in Table 1 as 
case 13, and were intended to duplicate the conditions of Fig. 9 of Ref. 102. 
The mass-flow result /-} = 0. 118 corresponds to 1. 577 x 10 lbm/sec. 
The thrust coefficient (with no correction due to ambient pressure) was cal- 
culated as 1 . 437 at X = H.O ( A /ft* = 34.10). This cor responds to a 

-3 

thrust level of 11. 12 x 10 lbf, and a specific impulse of 705 sec. The 
latter result agrees quite well with the line marked "theoretical 
performance" in Fig. 9 of Ref. 102. Apparently, this performance predic- 
tion is based on experimental data, and includes the effect of viscosity. 

Effects of Heat Transfer and Variable Wall Temperature 

Operation of a rocket device at low density implies that the thermal 
condition of the wall can exert a strong influence. Two cases were run in 
order to explore the effect of heat transfer on the flow. In the first of 
these (case 11) the wall temperature was equal to the reservoir tempera- 
ture all along the nozzle. In the second (case 12), the wall temperature was 
dropped from to T 0 / ^ between -4 and -2 radii, and was held con- 

stant thereafter: 

- r < * X - 4- 


— O. (c -f 0 . 4 - CPS £(x f 4) 

- o. Z , x ^ - z. 
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This latter case was an attempt to simulate the thermal environment that 
might be encountered when gas flows out of a hot reservoir and into a rela- 
tively cold channel. Cases 11 and 12 were the same as case 2 in all other 
respects. 

The wall-temperature and heat -transfer distributions, and the 
variations of the total -enthalpy flux (which are related by Eq. (3-31) ) are 
shown in Fig. 9, and typical profiles are given in Fig. 10. Profiles for the 
hot-wall case are practically identical to those of the adiabatic -wall case, 
but the cold-wall case shows marked differences. In particular, in the 
cold-wall case, about 20 percent of the enthalpy flux is lost in heat transfer 
to the wall. Also, the discharge coefficient is greater than one - -the nozzle 
can pass more mass flow than an isentropic, one -dimensional nozzle would, 
due to the elevated density level near the walls. The enthalpy profile in 
Fig. lOd has a slight bulge near the wall, typically of what happens in high 
Mach number flows due to viscous heating. 

The computer program described in Appendix A permits a more 
general wall -temperature variation of this type, where the end-points of the 
cosine variation and the lower temperature level can be specified as inputs: 

0-1 * < kTW/ 

'-'V , 

1+ TWZ | -Tw/ Z ( x-yTW l 

- +. CoJ rr ■ 

Z 2. \ xtWZ - xt i 

— TV 2 J X ^ XTWz 

A noticeable thinning of the boundary layer is apparent with the cold 
wall. Figure 11 shows the variation of the displacement thickness with 
distance along the nozzle, compared to the constant -wall -temperature case. 
The displacement thickness used here is defined by 


; VTWI 4 X < XTW2 
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Conditions on the axis can be calculated with an error of a few percent, by- 
using the isentropic, one -dimensional flow results at an effective area ratio 
found by subtracting the displacement thickness: 



Thus the present results provide a corroboration, without a priori assump- 
tions, that the displacement thickness is indeed the correct scale to use in 
making boundary-layer corrections to the core flow. 

Performance Characteristics 

The thrust coefficients obtained from these calculations are shown 
in Fig. 12, as a function of geometric area ratio. (No correction has been 
made for any effects that might accompany the transition from conditions at 
the given area ratio to zero-pressure ambient conditions. ) For comparison, 
the ideal value at the same area ratio is shown. These were calculated 
without any correction for divergence loss (see, for example, Ref. 103, 
p. 446): 
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The effect of increasing Reynolds number is to increase the thrust 
coefficient, as might be expected (cases 2 and 4). Increasing produced 

a slight increase in F (cases 4, 6, and 7), indicating that thrust losses 
due to divergence were not yet dominant in the range shown here. Decreas- 
ing the longitudinal radius of curvature of the throat (cases 1 and 2) and 
cooling the wall (cases 12, 1 and 2) both increased the thrust markedly. 

One interesting feature of the sharp-throat results is the rapid rise 
toward the maximum. For these configurations, the nozzle can be cut off 
at 5 radii with less than a five percent loss in thrust. These results illus- 
trate the importance, in these calculations, of establishing the correct 
initial state of the flow at the throat. 

In Fig. 12, the higher thrust coefficients are obtained at the expense 
of higher mass flows. The quantity F / is shown in Fig. 13. This 
quantity is proportional to the specific impulse: 


■m 


L 

4* fi 


If F 4^ is expressed in lbf/lbm/sec, and / 2- ' 
(for Y = 1.4): 


in ft/sec, this becomes 


[_sec£W& 5 ] = 0.0 ZZZ — • /IfT L FT_ /s£c] 

YV\ /) 


The results in Fig. 13 indicate that the specific impulse is a relatively weak 
function of Reynolds number and nozzle geometry, for the range investi- 
gated here. The strongest influence was that associated with strong cooling 
of the nozzle walls. 



One indication for design practice is that performance is increased 
by using a sharp throat, and by cutting the nozzle off at a rather modest 
exit -area ratio. This indication must be used with some caution, however, 
until more extensive calculations and further comparisons with experiment 
are made. 
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7. CONCLUDING REMARKS 


i 

i 


i The net result of the studies reported here is a computational method 

1 for finding a solution of the direct problem of low-density flow from given 

] reservoir conditions through a nozzle of given shape. Results found by this 

/;. method show good agreement with experiment. The equations used become 

identical with the thin-boundary-layer approximation at high Reynolds 
number, and can be used well down into the lower Reynolds number region 
| where the importance of the state of the flow at the throat makes the con- 

l ventional thin-boundary-layer model of limited value. 

j. Within the limits of the calculations made, it appears that sharp 

; throats lead to better performance, and that exit area ratios as low as 10 

j can be employed without serious loss in specific impulse. 

; .At sufficiently low Reynolds number and small divergence angles, 

the slender-channel equations do not have a solution in which the flow expands 
to supersonic conditions. Instead, the boundary layer completely fills the 
channel, and the solution takes on the character of a viscous, subsonic pipe 
flow. One implication of this result is that a relatively strong upstream 
effect of the exit -plane conditions is likely to be felt in these flows. Further 
clarification of this regime awaits experimental probing of the flow struc- 
ture and the application of more complex analyses which retain some of the 
effects discarded in the slender-channel approximation. 
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APPENDIX A 

COMPUTER PROGRAM DETAILS 


The program was coded in FORTRAN IV, and the G compiler was 
used. Execution was done on an IBM 360/65, Release 14. The amount of 
core used is approximately 45K bytes. A typical run may take between two 
and ten minutes, depending on how good the initial guesses at the mass flow 
are, how many axial stations are specified, and how much output is speci- 
fied. 

The required input data are described in comment cards in the list- 
ing (Appendix E). Most of these need no further explanation, but some of 
the restrictions on their magnitudes are given below: 

IBC: If this is read in as 2 , then XTW1 must be greater than XTW2 . 

IETAPR: The value 2 gives a convenient one -page output, with enough 

detail (51 values) to adequately define the profiles. 

Rl, THETA1, THETA2: The only nozzle geometry that can be handled 

by the present version is a converging cone joined to a diverging 
cone by a constant-radius -of -curvature section. The angles must 
be read in degrees, and should be restricted to values small enough 
for the slender-channel equations to apply. Angles up to 30 degrees 
have been used, but the magnitude of the error made is unknown 
at present. 

TW2: For cases where heat transfer is allowed (IBC = 1), the wall- 


temperature distribution must be specified. The only distribution 
allowed in the present program consists of two constant values 
connected by a cosine variation: 




The second level TW2 can be greater or less than one. 

B: This is a Reynolds number based on reservoir conditions and the 

throat radius. It can be expressed in the alternate forms: 


g = f>o r * _ / U W r* 

f*° J S '' 1 Jr o 


(A-l) 


Consistent units for the latter form are 




r° 

-r; 

<R • 

lbf/ft 2 

ft 

slug /ft -sec 

° R 

4. 97 x 10 4 ft-lbf/ slug-®R 

nt/m 2 

m 

kg/m-sec* 

°K 

i 

8. 32 x 10 m 2 /sec 2 -°K 


1 poise = 1 gm/cm-sec = 10 * kg/m-sec 

If the value of is less than about 600, it is probable that no solution 

with a supersonic core will be found. There is also an upper limit 
beyond which accuracy is lost, but its value is unknown at present. Runs 
have been made with E> as large as 3000, and it is probable that values 
at least to 10, 000 can be used. The loss in accuracy occurs for two 
reasons: first, as the boundary layer thins, it is described by fewer of 

the 101 equally spaced radial grid points. Secondly, the slip velocity at 
the wall decreases, and the integrals which are inversely proportional to 
this velocity are not calculated accurately by Simpson's rule over the 
equally spaced grid. 

Actually, for values of greater than 10, 000, a more conven- 

tional thin -boundary -layer program with no slip at the walls may be a 
preferable calculational method. 

XIPG: It is recommended that XIPG be greater than XMAX , so that 

mass conservation is enforced up to the saddle point. 
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ATEST: Present experience suggests that three significant figures are 

-4 - 3 

adequate; this means ATEST = 10 for V ^ 1. 3 , 10 for Y 1.4, 

-4 

and perhaps 3x10 for 1. 3 ^ V 1.4. If the mass flow is known 
exactly, in advance of the run, set ATEST = 1. 0, and AO = A1 . 

AO, Al; If the discharge coefficient is known, Eq. (3-27) can be used to 
make these initial guesses. If it is not, it is usually safe to assume 
C p = 0. 5 for a lower bound. The upper bound can be taken as =1.0 

for an adiabatic wall or a heated wall (TW2 > 1) and should be set at 
C_p = 1.2 for a cooled wall (TW2 < 1). The corresponding values of 
f\ can be quickly estimated from the table following Eq. (3-27). If the 
mass flow is known exactly, set A0 = Al , and ATEST = 1.0. Then 
the program will proceed immediately to the final solution, with no itera- 
tion on 

XCUT: An input is required here if the mass flow is known in advance, 

i. e. , if A0 = Al J set XCUT = 0 unless some previous runs have shown 
that it would be preferable to start the extrapolation through the saddle 
point at some other location. 

XPRINT: Subtract . 001 from the desired print interval -- i. e. , if profile 

output is desired at intervals of 0. 5 in ^ y. , read this number in as 
0. 499. 

DELX: The value 0. 1 has been found acceptable. It is not known how 

much larger this can be without loss of accuracy; smaller values can be 
used at the expense of longer calculation times. 

X0: A value of -5 is recommended, although smaller values can perhaps 

be used for B < 10^. 

XI , X2: If the throat is very sharp (i. e. , if R1 1), the transition from 

entrance to exit cones occurs over an axial distance equal to approxi- 
mately one throat radius. In such cases, it is well to set XI and X2 
equal to about -. 5 and +. 5, in order to follow some of the rapid changes 
that occur in this region. If R1 is greater than about 1. 0, set XI > X2 , 
and the small step-size region is omitted. 
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XMAX: This choice is dictated by the desired exit-area ratio. The present 

program has been run successfully at exit -area ratios up to 150. In 
some cases, an oscillation in may occur at X >10. Even 

when this does occur, however, the integ rated quantities (such as the 
thrust coefficient) appear to be relatively unaffected. 

ALPHU , ALPHT: Set these equal to one for perfect accommodation. 

XTW1 , XTW2: These are the stations that enclose the cosine variation in 

wall-to -stagnation temperature ratio. XTW1 must be greater than X0 
(preferably no less than X0+1 , so that the initial profiles have time to 
stabilize). XTW2 can be either greater or less than zero. No compu- 
tational difficulties were encountered when the difference XTW2 -XTW1 
was made as small as 1.0, with a temperature change A Qy - 0. 5 in 
that distance. 

The output begins with a printing of all the input values. Following 
this, the iterations on /) begin. For each value of A , for each y , 
and for each of the iterations on A? / jLy > a large number of quantities 
are printed. These quantities serve a dual purpose: some of the informa- 

tion they contain (such as the displacement -thickness variation) is often 
desired at more stations than those at which the profiles are printed, and, 
in addition, this output is useful as a diagnostic when the case fails to run 
to completion. Many users may wish to omit this output when making 
routine runs. A definition of all the program variables is given in Appendix 
F. This output is continued for each of the iterations on (\ , until A 

has been determined to the prescribed accuracy, ATEST . The final solu- 
tion is then calculated; on this pass, the output is augmented by a listing of 
the profiles at intervals of XPRINT . 

Some excerpts from the output of a sample case are given in Table 2. 
The first two pages of the output are shown, giving the input data, and the 
iterations used on the first two steps with the first A -value. Next is 
shown the page at which this first A -pass was terminated. The second 
A -value, 0. 095, was found to have too little mass flow, and the page of 
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output showing this determination is shown. The subsequent values of 
used were 

A = 0. 1025 
0. 10625 
0. 108125 
0. 1090625 
0. 10859375 

The last of these values is the one used for the final solution. Its first out- 
put page is shown, giving the initial profiles. The next page shown is that 
giving the iterations at X = -. 10 . The indicator IPG is set equal to 2, 
beginning with X = ^0 , since XCUT had been set equal to - . 1 0 on the 
/? -pass with = 0. 1090625. The next page shown is that beginning with 
X = 1.0. The saddle point is crossed between X = 1. 0 and 1. 1 , and 
six iterations are required to converge to the solution at X = 1. 2 . Finally, 
the profiles at X = 2 are shown. 

The quantities appearing in the columns of profile output are , 

l£ , O • Vs/ , \/ , , 4^ , M » and the five terms that enter the 

differential form of the streamtube-area relation, written in the form: 


VP A (A— 


/-N4‘ 
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rT 


xp 


ZY 


W 






DX 


PST 


^ (r-<) ^ 

— v 


2-^1 V ^ 
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IX 9 
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^ m 2 - ( y 
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CO^D'nJ 
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The last four terms should add up to a number equal to the first; these 
quantities are useful for studying the relative contributions that these dif- 
ferent effects make along the various streamtubes. 
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APPENDIX B 


DERIVATION OF THE GLOBAL CONSERVATION RELATIONS 
Continuity : 

Equation (3-15) is divided by P , and the result is integrated 


from zero to 

i ■ 

Jl^L + 

*"w 

0 

P 






? ( PIT 


V o 

The middle term is then expanded, noting that 


^ + \ ** (B-l) 


\a/ + ^ 

o p 




4L f* 1 Ai £^ L + p f* 1 i£- ^ j 

i 0 i o ) 


(B-2) 


+ zcrj f n dU£ ^ o 
J © 

o 

Division by j ^ Tl ^yi then gives the result shown in Eq. (3-25). 

o Q 

Momentum: 


In deriving this integral, it is simpler to replace P/q by 3 
Integrating the momentum equation from zero to one gives: 


f iff® = 


(B-3) 


/-I df , 

Zi <h- 1 B$l 


-±r M 


J 


^=1 
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The left side of this expression can be recast, as follows: note that 


4=^/^^^= f y* ^ + (' (Dl^) A(B-4) 

° O ° 

Use continuity to rewrite the last term: 

j= j'o<r|^ (»«)-— ^ (dA)! ^ 

Integrate by parts: 

U = | [a-yj 2 -^^ - j' 2. (vfu) | 

o 9 * ^ ^ ° 


-4 )[ trc 4 - r 


_^£ 


Expand these terms, and collect: 


4 - f ' y^-X)\ a — -<r ^ 6T ^ ^ | A? 

{ 9 * %y 'T'w 9v l 


(B-5) 


/ I 


- z 


^ w 


w 


J D(X Z y^ dy^ -b 


[Da-/6rC -vOL « i 




The last term here is zero, with or without slip at the wall. Using this 
result in Eq. (B-3) gives 
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(B-6) 
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4- % ci-x- 



2(X 
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Equation (3-29) is then obtained by multiplying through by 
and collecting terms. 


4Y«~w 

(y-0 


Energy : 


Note first that 


h = 

( B — 7) 

- j' n0|- J? ^ 

o o 

As before, the continuity equation is then used to rewrite the last term, 
and an integration by parts is applied to the result. This leads to 


_£ 

2X 


O O V 'W * 7, V. 


(B-8) 




*"v/ 


r 1 

2) 0 if y|_ ^ 


Thus the energy equation, integrated from zero to one, can be written as: 

A + - -y ^ l n* ^ = 


f g? ? 


^ ^ 


Sis' 




(B-9) 
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In order to obtain an expression for the total enthalpy flux, multiply the 
momentum equation by (J[ and integrate from zero to one: 




(B-10) 


Y~\ d-P 




J y ilX 


J 


f 6T — / 


J *\ 1 

L 9*[ / 


Note that 


_9_ 

2X 


j -a i = j ^ * 1 ! x h 


Again, the last term here is expanded, rewritten with the aid of the con 
tinuity equation, and the result is integrated by parts. This gives 



-*1 J)U 3 dr[ 



%-bai^Lt 

/ 



2!L + }yL 2J£ ) dvi 
) 


J ~£>LZ 3 *1 dr^ 


If this is now used in Eq. (B-10) above (multiplied by 2) and the result is 
added to Eq. (B-9), the result is Eq. (3-31). 
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APPENDIX C 


MATRIX - COEFFICIENT EXPRESSIONS 


The matrix coefficients appearing in Eq. (4-6) are: 
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APPENDIX D 


SERIES EXPANSION OF THE INITIAL CONDITIONS 

If the series expansions of (£ . 0 , P , and u/ are sub- 

stituted into the partial differential equations of motion, and coefficients of 
like powers of g^, are equated, the results are, for the momentum 
equation: 


~ < 3tt 3 ey + )' = ° 

41^ + +< a' j 

H + (v, a ■; )'= j 

& Y '^ 7r c ~ | ~ ^4 + 3^3 ) 

(D-l) 


4- U/^ 6T 4 4 4 lv^ LX £ j 

7*-, + fx«V)'+ (^OT+uL)' = 

= f -rj t7lf,ti 4 ) - <rJ-r } -2iZ? +< iZ r ' 

4 \a / 3 ^ + w 4 x 5 ' + u/ r -t- Vv/j ! 
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• 


For the energy equation, only the leading term has been worked out: 


T- < s n ^ )' + 2^ («•' y 


(D-2) 


The continuity equation is used to find l^/ directly from (JC , (£) , and 


P : 


1 Vv/„ = /a/- 2 ) J y * 2,3,4. 


(D-3) 


' *1 

j +3^' 7T [ UC^^d-x 


As noted in Section 5, the coefficients 7T, are found by enforcing the 
global conservation of mass. 

Results for the first few terms are given in Section 5; the remain- 
ing terms that have been worked out are 

e 1 * 3 1 tjf ^ - T ("'*) ("’Oj 

■t-<i^'/l 2 6C j- j -zaf,-^) +zs(i->i 4 )- 8 (t- z^ z ) 


W ?- ^ 3 - f n r - £ V 


5-400 


+ K/) v /4 Z gC -|.^4/C?^ 3 - Z^ 1 J + 32 < (^-X 3 ) 


TT r = 


/j3d 

27 OO 


rj s t 3 + + 22 22 


5* rr / 
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*'<*.) - (rtf** 1 + -g? o-*6 - £ M'J 


+ f 0-0 -~r ('-\' 1 ) ■<- 


8 loo 


LLtSo ^ ^ ) 


+• 


A 3 8 Z K U , |- y (i-^) + ~^ (t-\ 4 ) ~ y 2 M 4 )* y ^Vyr 0'^°)J 


+ [z 4 o (,-rf) - IL 8 (<-\*) + (,-y[£)J 


+ 


[- 


/Z <5 I 30 /& 


3 (r-i) zj 8 \i8t>os 


( r„')Ve s + to'f w* 


< 4 r 


- /44 %l A Z 3 K* - 7St> A j 


- - / 3 ._ (tt'^ ai 


(rtf A 3 3 “ K* + Z 4 /? 1 5 J^t 2 + <28 A X { 


US’ 


3 


-i- V^L crj £ tt; = 
0 ZY ^ 


sir 4 


3 (fr-i) g-JB 7 ' 44 z<? 


fc)Vs J 


+ ~3LtJ -C4fi K? 
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APPENDIX E 


PROGRAM LISTING AND FLOW CHART 


■; ciaa . . _. ; _. i 

KUL tAVIJ] 2 

CCMMLN A.AbULPi . AB01-P2. AaCPGl , AtIL)fG2. AEFF, ALPHT , ALP FU ,A MAC H , ._ .. . 3 

: *AMACH2,AMM2,AMM2A , AMM2A V, AMM28, AHM2C, AMM2MN , A MM2 l , AM S C . AM2LS T , 4 

! *AB£i-l»AItST , AO , A1 , Al 1, A l2 t A21, A22 . .. _ . 5 

; COMMON d, BOTTOM. aKb,BKl,3RU,efll2,fcR2,BK2i,BK22, dll, U12,d2i,B22 6 

; CCMMLN C, CFKGT, CUNOI.CPKU ,CU,C 12 ,C2l ,C22 _ 7 

1 CCMMCN G,L)ELETA,UELA,UtLXS, CfcNGM , CF BKC A , DL P 1, OL P2« OL S TR , DPI) X , DP G 1 , 8 

»UPG2 .L)PLAST .GPNM1 .JPNuH . 0 FNG MA .UPNO M B . OT ON . CUC N , CU ONS C , DVUN, PI, 02, . 9 

! *L27LN2» C2G0N2 , U A Y ^ 10 

j CCMMCN EMll tEMW i EM21 i EM22 ,EIJ Eu2 .EilMAX 11 

j CCMMUN FACI.FNl.FN2 12 

l CCMMCN GAMMA, GAMMAP.GfcNRI , C l ,G 1 C , C 1 i , L 12, G_l 2 xG l A ,G l 5 ,0 l 6 ,G 1 7,GL8, 13 

* G 1 4 ,C2 , u2 0 , G2 1 , G3 , G4 , G4 BGM i , G5 » Go . G 7 » G8 » G9 14 

C OMMON FUEL X.HF L UX. HIK . 18 

, CCMMCN 1 fcc, IETAPk, IGuTA.lPG, ISw » l T tK, ITEkA, l X , I X ST AR , I XTHW 16 

COMMON K ,K1 ,KPl , NX ,KX0, KX l 17 

i COMMON OC.MUKb.NPKO.NKUN 18 

■ COMMON CF;MPFt, OMEGA. Cl, G 2 , C 3 t£_4 , 0_5 ... 19 

COMMON P ,HAK I l ,PAKT2 ,PUA, FG , FGNP1 ,PGKI ,PG1 ,PG2 ,PFW , PI3 , PNt« , PN1, 20 

•j ... *PN2 , PKj.PS.Cfi.2t PS 1 , PX TK AP _ ... __ . 21 

j COMMON w , GO ll,Cwl2, CC2 1 ,C022,Cll,C12,021,C22 22 

i CUMKCN KA0UtGiKIPa2,KT2827 l KI26M.l..HA»HliO»Mli.U i KU .23 

i COMMON SAMb.SA 1,SC1,SGDX, SFER I , SICeAR.S IGUP.S IGMA*,Sw2 ,St.i, SW4 24 

i COMMON. J A NtlAK, 1ANTHP,TAU,TCF,TCFTCP»T FCQ££ xTHEI AP , TUtT Am, T Hfc T A 1 , 25 

; *IHETA2.TCMF ,TCiP,TltKM, Tv»2 26 

| jCULMCA GCF,oOFTGP,oTEKM f U2l rAT.y/ ... _ ... 27 

•] COMMON X,XC0T,XHN,X1PG,XMAx,XPLASI,aPHINT ,XJWl,Xl»2,XC,Xl,X2 28 

COMMON Yl,Y«tY3,Y4,Y3,Ye t YT.VE.YS 29 

COMMON 21, Z10, 211, 212,213 ,Zl3,Zlo,Z2, 223, 231,23d, 232, 234, 233, 24 1. 30 

*242»24J,Z5,27,28 31 

COMMON ET AllOl), £111 10 l ), t 1 2 ( 1C 11 ,E21( IC1) ,t22( IC1) ,F ItlCl'l , 32 

.. *fe2l.lCl) ,*311011 ,S4l.l011 ,MUOl) 33 

COMMON II 101,21 ,01 101,21 34 

C 1 Mt NS l UN AKCFi ( 1 0 1 ) , PGF.O I 1C 1 ) , CUNC ( 1C1), SHERI lUll,GtNK(101) 35 

; ECO 1 V AOE NoE ( ARCH , t2 l ), IP OH U , E22 »,( CCNU ,F l J ,( SH tk , F2 ), 1 Gb NK, S 4 1 36 

COMMON/ERRO 1M/ CHECKS14) 37 

! I N 1 ECLR*2 CHECKS 38 

CHtCKSm =1 ... 30 

CFECKS(2)=1 AO 

CHECK S( 3 >=-10 41 

i CHECKS! 41 =1 42 

.! CAOO LVJCHK 43 

500 CALL LO E Ak 1 A ,U I 1 01 , 2 ) I 44 

3 CALL 1NPC.I. _ .. . _ . ... 45 

1 CALL O', IT 46 

47 

dtGINMNG OF A ITERATION LOOP 48 

C 49 

1 CALL START 30 


: ,X 


•v 
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: I I , jog 




BEGINNING UF X I NT EGKAT 1 CN LOUP 


PGKT = 1 .U 

IPG =4 

IK X.GI • X 1 PG ) IPG=1 
4X1=0 
I XS1AK=1 
iSh = l 

5 UC 4 L= 1 , 101 
UL, 11=0(1., 2) 

4 I ( L *1 1= I t L * 2 ) 

IT£«=l 

PG 1=0 FU x*PGk T 

P0*=Li.Ol * PC* 1 

GC TO (301,3051, IXST Ak 
1.0.1 IHK1L^1»5J DQ. ID 105 

IF (CPCX .GT .UPLAST J 1XSTAR=2 

21 IF(X.LT.Xl) GO 10 23 

14*2 .... 

Kl=Kl*10 

01 LX =Dc.L X / 1 0 , 0 
24l=234*CELX 
PGKT=PGKI**0.2 
GC IC 23 

22 IF (FCCIKl.lO) . NE.O) G O 10.23 
IF1X.LT.X2) GO TO 23 

14 = 3 . 

KI=KI/10 
UELX = Ufcl-XS 
24l=234*CELX 

.. RGRJ f£Gr< 1**5 

PC l=l)PO X*PGK T 
PG2K.01 *PC1 

23 K I = K I + 1 
CPLAST = OPl)X 
AF2LST=AFP2 

UPCX=PG1 

X= XL ♦ PLUA 1 IK I ) *UELX 
KX=KX+1 



PEGINMNG OF P ITERATION LOOP 


MCLlX = 0,S *Dfe LX 

Xh*=X-HCELX 

GO TC (41,42,14), IXTH4 

41 IF(X.LE.XTwl) GO TU 14 
I -XT PJX-2 

42 IFIX.GT.XT42) GO TO 43 


51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
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3 

i 


j 

j 


• 3 
•i 


• 1 


UltIAi«=ia*Ui*CQS(.G5*lX-XIJ»JJJ 

GL TC 14 


101 

102 

103 

4 


rHETAh=Tw2 


104 

j 


14 CALL G£CM (.XHJtl : _ ... 


IC5 



SIGEAP»SlGf4Ato 


106 

■ \ 


] ANdAR* IAMHP 


107 

' 


PHR=P*JPUX*HIJELX 


1C8 

•i 



LLLliiziiU ... ... 


109 

' i 


fc 12 11 1*0 


110 

• - A 


t2111)*C 


111 



1:22(11=1.0 


l 12 



FI ( 1 ) *0 


113 

j " .j 


F2( 1 J=0 


1 14 



YL=235*51GMAW**2 


115 

~\n;| 


Y2=2.0 *Y1 


116 

' •- -i 


Y3=4.0 *Yl 


117 

. 1 


Y4=F.C *Y1 


1 18 



Y 5 = 2 2 4* S I OMAR 


119 

. A 
■■ ' -j 


Y 6*2 4 1*4 1 GM A» 


120 

" j 



Y7=YL*PR . 


121 

< 


Yfl = 4.C * Y 7 


122 



Y 9* 2 7* CP LX 


123 

} 


GC o 1*2.100 


124 

J 


. Gl=PHta*fcJAtU .... ...... . . ... 


125 

4 


C2 = 01 *m (L 1 


126 

■ J 


GJ=Y5*T(L.1J 


127 



G4=I ( Li l) **232 


128 

„ _ j 


G5 = 233*r < Li l 1 


129 



Go = CMtG A* ET AIL) 


l 30 



G7 = T4L + l.l)-T(C-lf 1) 


131 



G8=G6*G7 


132 

•*. 


G4=4.C *t FAIL) j ........ 


133 



G1U=C9*T<L, 1) 


134 

i 


G11 = U(L*1 .1 |-Ul-l,l) 


135 



G 1 2= I1L , l)**0MfcGA 


1 3o 



G13=£TA(L)*G12 


137 



G14 = CEL£IA*UL,1) 


1 38 



G l 5 = 2 . G *E 1 A ( L I 


139 

* ’’ l 


Glo*G 15*7 l L. 1 1 


140 

■A t \ 


G 1 7=UtL 4* I ( L . 1 1 


141 

• ' ■•! 


Gia*Gl*U<L.U 


142 

i 


G19=GELX*W(L) 


143 



G2C = Yt*I( L, l) 


144 

■ • l 


G21=»9«tIA(L) 


145 

' /v • ' 

c 



146 

. •! 
■ t 

c 

CALCLLATIUN UF THc MATklX lLEMEMS. 


147 


c 



143 



A l 1=-G2/G3+G4*< G5 *08*010 2Y4 


149 

■ ■ j 


A 12 = G4*G4*G U/Y4 


150 


■A 


; r'A 

r U 
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I ■ 

h; 



it>i 

152 

153 

154 

155 

156 

157 
15B 


£Xl>C2/iL}-C4*i(i^t£S=iaO)/.Y4 _ 159 

C 1 2 =— A 1 2 160 

£21=-A2l 161 

C22=G2/G3-G4*1G14+G8-Gl6l/Y6 162 

C 163 

C1=G1*( Y5*UIL, l ) *♦2-019*0 1 1 )/G2C-G21/2£ 164 

♦ ♦G12 * 1DELE1A*G1 1 + G 1 5 ♦ 1U 1L»1. 1 1-2.000* 0 IL « 1 1 *U 1 l- 1 . 1 ) 1 ) / Y 3 165 

C2=Gl»i Y5*ll(L, 1) *T(L , 1 1 -G 19*67 1 /G2C+C21 *0 1 L . I I /GAM PA 166 

* ♦Cl2*lD£LtTA*G7+Gl5*(T(L*l«l>-2^CU0*Tit«JJ*riL-l,l) I ) /YB 167 

CCU= 622-C2l*E12lL-ll-C22*fc22IL-l) 168 

uC12.=-B12 + Cll*tl2lL-l l+Cli*E.22(L- U 169 

tC21=-821+C2l*fcl lll-l) +C22*fc21l L- 1 1 170 

Bll-C U*i-l 11 1- 1 1 -Gl Lg.S2Uk-J.l_ _ m 

CEXCM=CCll*Ca22-Q02l*g012 172 

ttll=«iill/DENQM 173 

0 1 2=GQ 12/0EN0M 174 

w21=CiZl/BENQH 175 

C22 = CC22/CEMJM 176 

pakti =cii *f hi-h+ci 2 *f 21 i-i>+ui ...... . 177 

FART2=C21*F IIL-1 1 *C22*F2 IL- 1 )*L2 17B 

£11 111 1* All* 612.*A21 ... 179 

E12(U = Qll*A12*y 12*A22 180 

£2ilU = C21*AUtC22*A21 181 

E221Ll=Q2l*A12+022*A22 182 

FltLl ?1L1 1» PART1 +412*PA8T2 183 

t F2(LJ =(.2l*PAKTl + G22*PARI2 184 

CAtL 0EUM1 X 1 ... 185 

S*«2=5 1GMA***2 186 

TCJP= C*7 SIT HET A P 1 * T 1 1 0.1 * l 1 •♦CHMPH 167 

d(JITOM=CFdOT*l P+UPL»X*OELxi*SlGF.AH 188 

LCF = T CP* UCFT GP/BCTIC H _ 1 89 

TCF = ICP*U.F TUP/8UI7UM 190 

FN 1 = U( 101 .JJ-U11CC . j.lJ-F.H 1G.C ) .... . . ._ 191 

EK l l = tl 1 1 100 1-12.0 ♦UCFl/UCF 192 

ePU = tU(l.QC 1 _ _ 193 

lFllEC.Nfc.il Gu 10 7 194 

FN2=— 12.0 *1 HETAK/ TCF 1+T 1101 .11-T11C0 .11-F211CC 1 195 

EM2 1=E2 14 100 1 196 

£F22=E221 100 1-12.0 ♦ICFI/TCF _ 197 

GU TU 8 198 

7 FN2 = IU£i ,11-JIiOC » 1) -F2< ICC 1 _ 199 

EM2 i = E2 1 1 100 1+4.0 *PR *U 1 101 , l) /UCF 2 00 


A21=G13*Gll/Y2 

A22=-G2/G3+G4»1G14+G8+G16 )/Yh 


81l=G18/017+G13/Yl 

_ £12=C 

82l*C 

ti22=£ia/C17tG13/Y7 
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8 0ENCM=EHll*EH22-£M2l*EMl2 
- UtiCl .2 I = 1EH22»EN1- EH 12»FN2)/UCM> 


TU01 ,2 » = ( EMI l*FN2-£M2 l*KNll /UENOM 

U-rLlijLUZJ . 

T l=T ( 101.2) 

Eiitiai)?ui/Ji 

fc 12( IUI)=UI*E1UICI) 

E22(1C1)=T1/U1»*2 

CO S 1=1.100 

u< ici-i . 2 ) =e in iui-n *ui»£i 2 ( ioi-i i ♦ ti *ei ( l oi- n 

Ul=G( 101-1,2) 


IF (U .LE.O.O.ANC.X.Gl.O .0 ) GC TC £$_ 

K 10 1- 1 ,2) =E21( 101- 1 ) *U !♦£ 22 ( IG 1— I ) *T 1 *F2 ( 1 01 - 1 ) 


tllIlCl-l)=UI*ETA(l01— I)/II 
E22( I0l-I) = EJAU01-I )*TI/L1**2 
E12(ICI-I)=UI*E11(101-I ) 
Eim2r*iJi .*t22(2l . .... 

SA ( l ) = 0 


I ♦El 2 (2 1 ) ♦ A j__C 


SA( 2) =2 23*1 a.O *E 1 1 ( 2 ) — E 1 1 1 3) I 

02fcl AT=A.C- *E12 ( 2) 

00 10 l=3 .1 01 
SCALE=A— 2*MCU1L«2) 

IFIL.bg. 101 ) SCALE = l 

ETTHU2 = EI TB02*SCA Lfc*E2 2 ( L I _ 

02EIAT=U2ETAT*$LALE*E12(L) 

HELliX = HFLUX*U(Lj2) *( ETAIL lt£12(L I )*SCALE 
10 SA(L)=SA< L— 2 ) > Z 1 3* ( fc 1 1 ( L— 2 ) ♦ A • 0 * E 1 l I L- 1 ) ♦ 1 1 l ( L > I 
AMSC=ZA3*U2ETAT 
AMM2=C.5 *Z7*Z13*fcTTBU2 

AMM 2AV=0.S * ( AMM2*AM2L ST ) 

IF ( IGCf A. EG. 0. ANU.AMM2.LT .0.51 GO TO 2C2 
PNE*=A/< SA( 101 )*SW2) 

PXT«AP=P*UPOX*UELX 
AKCF I =0 
SHE HI =0 

..... GUM)l = Q . 

GENHI=0 

PSG*2 = PNfc«*Sl»2 

SGCX=SIGMAw/CELX 

HFLUX=PSGw2*Z13*HFLUX 

THC CEF = P SGh2* ( 1.0 ♦Zd*AMSC) 

0TE_RP=(Uil0L»2)-JJ( 100. 2 ) + L( 161, l)- 0( ICC, 1) ) 

TTERM=(0.5*( Tllol ,2)+T HOI,!)) ) **LMEGA 
0= ( T T ERM*( (1(101, 2J-TC 100,2)±T(1C1,1)-_J( 100,1 ) ) /PR 
* ♦GIERM* (01101, 2)*ui ICl.nh )/(6*Z33) 

.. E1K1C1 ) = PSG*2*.SA(l0l . I 
E 12 (101 )=0 


201 
2 02 
203 
2 OA 
2C5 
2 06 
207 
2C8 

209 

210 
211 
212 
213 
21A 

215 

216 

217 

218 

219 

220 
221 
222 
223 
22A 

225 

226 

227 

228 
2 29 
2 30 

231 

232 

233 
23A 
235 
2 36 
237 
2 38 
239 
2 AO 
2A1 
2A2 
2 A3 
2 AA 
2A5 
2A6 
2A7 
2A8 
2A9 
250 
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. t Li 1. 1.1 = 0 251 

UC 2 L = 2 » 100 252 

ElltLl =PSGk2»S5{L) 253 

2 £12(L)=0.5*(SGDX*(T(L,2I • S4 ( L WE I A 1 L ) ) * ALCC I S 3 (L) /E 1 1 ( L 1 1 ♦ W I L 1 1 2 5* 

GAHiEAP=liAMtlA*j?Hlt .. . . . .. 2.55 

EACI = i!8/(27*B*SIG8AK) 255 

SJ.GO£=S1GMAj<*OPUX 257 

AKCE(11=0 258 

CONDI 1) =0 _ 259 

SEEK ( 1 )=0 260 

GfcAPJLL=8. 261 

PGR 01 1) =6 262 

CC 3 L=2, 100 263 

CUCN»lUl LU ,2 1 -GlL-1.2) *U l L ♦ 1 > 1 ) -U IL— 1,111/2 34 269 

UVChi* IE 12 1L ♦ 1 ) — E 12 ( L- 1 1 ) /2 33ETANBAK* 10. 5 »(U(L,2) ♦U(L,1») 265 

» * £T A I L i »DUDN ) 266 

0 TOCsS.l UL*1 .2) — TtL.-l »21 ♦ ! (I ♦ l , l l-T IL- 1 , 1 1 1/ 234 267 

020UN 2- ( U I L + 1 , 2 I -2.0 *L(L,2I ♦U(L-1.2) , 268 

.... * - .tlUl.ti. 11-2.0 *UU, 1>ECIL-I,.lt 1/242 269 

D27UN2=(T<L* 1,2) -2.0 *I(L,21 «T(L-1,21 270 

* ♦Kl«l. 11 — 2.0 MIL. HeUL-tI. 111/242 271 

APAO = 0.5 ♦KT2GM1*(U( 1,21 / SCPIITI L,2) 1 272 

* ♦ U ( L, 11/ SORT I T ( L, 1111 273 

AKA(.E2= AMACH**2 274 

SAC£-=ln.j2J ttIMJLl*C;MfeG/*C7CN 275 

V= E 12 IL 1 ♦ T AN T HP *E T A ( L 1 *U( 1,21 2 76 

AKCE(L1=CAMHAP» (V*(1.C -E T A 1 L 1 *DUUN/U1 L , 2 1 1 *E 7A1 L 1 *U VDM 2? 7 

* / L ( L ,2 ) 278 

IGH7 = EACJ *1 lii 2 1 ..**232 . . .. 2 79 

SHfcKlL»=-TOME*(OUON*SAME+HL,2) * E I A ( L 1*D 2LDN2 1 / AMACH2 280 

CUM! 1 LI = 1 TOME * 1 0 TON* SAM t* 1 1 L ,2 1 *E I A JL | *D2T DN2 1 *U ( L , 2 1 / 281 

* (1 1 L , 2 1 *PK11/AMACH2 282 

_ r,ENK 1 L1 =2 .0 *TCMF »U(L .2 1 ♦tTAd 1 »OUCN » *2 / AM A CH2 _ _ 283 

PGR0m=(l.U — AMACH2 1*ETA(L)*S1GCP / A EACH 2 284 

SCALE=4-2*MGU(L ,2 1 _ 285 

AK CE 1 =ARCH l ♦ SC AL E * ARCH ( L 1 286 

SEEK 1 = SHE Hi ♦ SC ALE *SHER (El _ 287 

CCNCl = CCNunSCALt*LUNI)( LI 288 

3 GENRl=OENRl ESCALt*GE NR (Ll . _ 289 

ARCEC101 1=0 290 

APACh=0.5 *RT2G?U*(U< K1 ,2J / JUKTJ T.( i_CJL ,2j ( 251 

* ♦ UllCl , 11/ SwPimiOl , 1)11 292 

AMACE.2=AMACH**2. _ 293 

PGRLll 10l)=S 1G0P* ( I.U-AMACE2 1/AMACE2 294 

ARCHl = GAMPAP »TAN6A R 295 

SHE R ( 10 1 1= -FACT*B*Sk2*l27*OPUX/2«EPHk*U 1C1 ,11* 296 

* (UUOl ,21 -UllOl ,111/17(101 ? I 1 *UEL X 1 1 / AM ACH2 257 

CUCNSU= 11 U( 101 ,11-UUCC , It E L ( 10 1 ,21 -UllOC ,21 1/233 1**2 298 

CCNC1101 1= FACI*(UJiOl ,2i ./T1JC1 ,21 l*l_-2.0 *7(101 ,11 299 

* **UMEGA*DUUNSQEb*SR2*<-27*U(10l , 1 1 •GPGX/ GAMM AePHW*UI 101 ,11 3C0 
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* *U(1UI .21 -T1101 .Lll/ milll .1)«UELX)1)/AMALH2 

GfcNk ( 10 1 ) = FAC1 *U 1 1 0 1 i 2) *2.0 <1(111 .21 **Z32*DUDNSG/AMACH2 

ShfcB l«Z13*lSHERI*SHEIM 101 )) 

CCNC1 = Z13*ICCN0H-C0NUI 101 )) 

-CENK 1= 2 1 3 *1 GENK 1 »GENR ( 101 )) . _ 

OPMjM=ARCMl-SHER I-CGNU l— G£NR 1 

PGN F l=QPNOM/(SIGU Ak* ( AMM 2 A V— C .3 1>_ 

0 THE FOLLOWING STATEMENT It USED TC SUPPRESS OSCILLATIONS IN DPOX 

. 1 P ( X . GT .10.01 PGNP 1= (P GNPUCPLA ST ]*C.5 

AEFF=2.C * A* T ( 1 . 2 ) / ( PNE h*L I l ,2 ) > 

CLS I H=S 1GN Asr. SURT1 AEFF } _ 

IF(MCC(LC.48).EC.O) MR I T E (6 , 2 17 1 A 

Jrf 1 Til 6.2231 _ X.lPlX.PSNPl.PXTKAP.PNEw.on ,2 ) ,T{1,2), 

* AkCHI .SPERI.CONUI .GENR i .DPNOM.AMS0. 

* hFL UX. S IG8AR.SlCM Aw.TANBAk . I ANT HP . AMM2 » 

* 0(101 til . T ( 1 ( 1 ,2) . AEFF > CLS TK tC 

LC=J.<L<4 .. ... 

GG TC 151. 600. SL. 620) . IPC 
__.6HC GO TU 1.68_L».6Q2.i 19 1 , ITER , . 

6C1 AMM2A=AMM2 

OPNLM A =UP NUM .. 

OPG 1 = PGNP 1— PG 1 
I T EK = 2 

0 PU X* FG2 

.... GC TC 14 

602 AMM2U=AMM2 

flfhUfB=DP<NUM 

ITEF-3 

lFlISts.EQ.21 .££_ T 0 599 

6R 1= ( OPNUMb— OPNUM A ) /(AMM2B— AMM2A)+DPNMl/ (C.5 -AMM211 

8R2=CPNM1-PPNUMA*0PNM l*AMM21/<0.5 —A PM 2 1 ) ♦ ( OPNLMb— DPNUMA ) * AMM2 A/ ( 
♦AM M2 8— A MM2 A) 

AMM2L »BR22BR1 .... 

UPCX= (PGA*( AMM2B-AMM2C ) + PC-2 < 1 AMM2C-AMM2A ) ) / ( AMM2B- AMM2A ) 

GC TC 14 

620 IF ( ITER.GT.2) GO TO 62tt 
IF (I TER.EQ.2) GO TU 625 
ITER=2 

PNl= PN£k 

CLP1 = FM-(P*PG1 *CELX) 

0PQX=PG2 
GC TC 14 

625 FN2=FNEW 

1 I Ek= 3 

. OkF 2= PN 2- ( P.+P.G2 *O.EL X ) 

626 PG=( (PGl*PN2-PG2*PNU-P*(PGl-PG2l 1/ 

* ( (PN2-PNU*0ELX*(PG1-PG2) ) 

IF ( PG .GT .1.0 I PG=l.0E-02*PG 
AtsCLPl = ABS(OLPl) 

ABCLP 2= AbS(0LP2) 


SCI 

302 

303 
3C4 
305 
3C6 

307 

308 
3 09 
3 10 
311 
3 12 
3 13 
314 
3 15 

316 

317 

318 

319 
3 20 

321 

322 

323 

324 

325 

326 

327 
3 28 
329 

2 30 
331 

3 32 
333 
3 34 
335 
3 36 

337 

338 
3 39 

340 

341 

342 

343 

344 

345 

346 

347 

348 

349 
3 50 
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if lAMINli AbULPi.AtJDLPZi .EL. ABOLPll GO T(j 627 351 

CLP 1=CLP2 3W 

EGJjlP.G2 _ 353 

PN1=PN2 354 

- . 622 .££2 ?PG. . .. 355 

OPDX* PG2 356 

m JO 1.4 ...... 357 

626 P62 = FNEw 358 

ULP2.»PN 2rrXP »P ..G.2*ll£LX 1 . . . .. . 359 

Iff ABS(OLP2)/P.LE. l.Ot-O'l GC ID 625 360 

I T EK= IT EB+1 361 

If l 1 ItR.LE. 10) bU TO 626 362 

*&J.Uiiu2m. .... 363 

If (1GCT Aj EQ.UI GO TO 304 364 

GJJ_L0_1A _ 365 

625 IMX.GT.xIpG 1 IP G= 1 366 

GO IC 19 .367 

51 IF ( 1 1 Ert .0 T . 2 ) GO TO 1 8 368 

IPLUtR.£a.2J GJ TO 15 . ._ 36-7 

l TER = 2 370 

U£G.l=Pu£lPl-PGl _ ...._ 3 71 

CPCX-FG2 37 2 

GO 10 14 273 

15 UP G 2 = PG NP 1- PG 2 374 

II£R=i ............ 375 

16 PG= ( PG 1*0 PG 2— P62*U Pul I / (DFG2-CPC11 3 76 

If IPG.GT-l.O 1 P G= 1 . OE-O 3*PG 3 77 

A ELPC 1= ABSIOPGU 378 

AB0fG2r .AQ51UPJG21 . 3 79 

IFIAPINIIAHDPGI t AB0PG2 l.fcw.ABUPGl) GO TC 17 380 

UPG1=CPG2 381 

PG 1 = PG2 382 

17 PG?=PG 383 

UP0>=PG2 384 

GO 10 14 _ 385 

18 CPC2=FGNPl-PG2 386 

lf( A8S({)PG2/_DPLAST)._Lfc.l .OE-03) GC TC 19 387 

ITER* I TER ♦ l 388 

LF.J lTfc»«.U. ULLJit„UL L&„ .._ 3 89 

«k 1 IE ( 6 i 2 10 1 390 

1£1 1GCTA.EQ.01 jGO 10 304 . . ...... 391 

GC TC 13 352 

LS.P.-PJtTKAP . ... . 393 

hRIIE(6#224) 394 

K6B.l»teu*/im,A£I _ 395 

1HAPP2.LT.AMM2MN) AMM2PN-APP2 396 

CO 11 L=l,10l _ ■ . 397 

Ml»=E12(U 358 

ll S3( L )-E IK L 1 _ 399 

If ( IGOTAiNE.Ol GO TU 90 400 
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... - ihxi:£XJi«fir^a.Qj GQ.ja.2i 401 

GO JO (9i,92), IXSTAR 4C2 

21 IF t P «GT .0.1) QU IG ? . ■ ... 403 

GO TO 304 404 

.22 _i.F_LU_£UXjGJ .oPLAST) ,GCL T.C .3 _ . . _ ... 405 

IFiKX-XXl .Nt.U GO TO 303 40o 

KX J=iK X 407 

iF(KX^KXC.GE.3) GO JO 300 4CB 

- G0_IC_5 .... . .. _ . 409 

303 KX 1 -KX 410 

KXC=KX 411 

xai = * 412 

GO rc 5 413 

3 C 0 bfllIECO,203) 414 

GQ LG 3Q4. .. .... 4 15 

302 *RlTElo,204) AMM2 416 

304 A 1 - A ... 4 17 

XCLJ=X-l.Ol *UeLX 418 

GL TC 23 . 419 

27 AC -A 420 

#Kirfc(6,20/l _ 421 

26 IF ( AES UO-Al ) .L E.ATEST I IG0TA=1 *.22 

GO TO l 423 

424 

END CF A ITERATION LOOP 425 

426 

13 if.IMOK.tt-Nt.Ql GO TO 500. 427 

STOP 42b 

SC if (OPDX.GT.O.O) GG 10 9S 429 

F.TR = hTR4i«»DeLX 430 

IF (NFPa.EC.O) CALL UOTPuT 431 

GO JO ( 23,94, S8. S3) , IPG 432 

23 LfiX.Gt.XMAX) GO TO 13. .._ _ .. ...... . 433 

IF(X.GE.XCUT) GC TO 931 434 

IF ( < AMM2-AMM2MN) .GT.0. O.ANO. X.GI.-0.5 ) GO TO S31 435 

GC TC 5 . 436 

231 l PG = 2 43/ 

AMP21=AMK2 438 

CPNF l=UPNOM 439 

GC TC 5 440 

94 IF(AMM2 .LI .0.42992. ANO. DPNLM.GI .0.0 ) GC TC 95 4Sl 

If IX.GE.XKAX) GC TU 13 442 

GC TC 5 443 

95 I SM=2 444 

._ GC JO .5 ... .. .. 445 

592 I PC = 3 446 

EH 11 = ( OPNUM 6— GPNUM A ) / < PG2-PG l ) 447 

bfil2=-8Rll*PGl*0PNUMA 448 

BK<l=iAMM2B-AMM2A) /(PG2-PGI) 449 

BH22=-bR21*PGl+AMM2A 450 


6 K6 = S I GO A K*(AM2L ST+8R 22— 1. C (-2.C *BKll 

DPUX=-(8ke* SuRT ( 8kd**2 +8 .0 »8K2i*bK 12*SIGBAR ) ) /( 2.0 *B R2 l * SI G B AR ) 

£li2=iU?0.X. . _ _ ..... . 

GC TC 14 

S3: If.U.LJ.XMAXJ GG TO 5 . . 

GO TO 13 
.99 «R JJ£ 16, 2 C5 ) 

GC TO 13 

tt. J.= 1.L1-J . . . _ 

MR 1 TE ( 6 ,206 ) J,0(J,2) 

GC TC 13 . . 

203 fURMATi • JPOX IS UECREASING fCK TCC MANY STEPS' I 

2.04 fO BMA Ii' AH M2= » . 1PE14.6) . . .. 

2 C 5 f CHMAI ( ' OPUX.GT.O.O ON PBLfllE PASS') 

2 06 FORMA T 1 ' NEGAT IVE VELOCITY ENCOLNTEREO AT G 1 ' , I 3 . ' , K+ 1 1 = » .1PE13.6) 
207 FGPMATC TOO LITTLE MASS FLO*') 

230 FGtiMALLI.' JQC MANX LTEiJAIICNS CN P’ I 

217 FORMAT i • 1A=', IPE16.6//J 

2X3 F CRfLAliJ)££i2t 2j4AiiP6£.lJj ,6«3J/loX,6El6.6) ) 

224 FORMAT! IX) 

..fiHC . .... ... 


451 
4 52 

<*53 

4 54 
455 
4 56 

457 

458 
4 59 
4 60 

461 

462 

463 

464 

465 

466 

467 
463 
4 69 
4 70 
471 
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cccc 


4 72 

SUBROUTINE GE0MIY1 <.73 

_J4£At O.AYI2J . ... . . 474 

CCPPCN A,A8ULP1, A60LP2, ABCPG 1 , A C0PG2 , AEFF , ALPHT , AL PHU , AMACH , <.75 

*4 PAC HZ , APPZ,AMM2A, AMH2 A V . APPZU, AMP2C, AMM2 MN , AMM2 1 , A PS C, AM_2L S T , <.76 

♦ ARCH 1 ,A TE SI , AO. A1 .All , A 12 ,A21 ,A22 4 77 

CCPPCN B,80TICM,BRB,BRl,8Rll,8R12,BR2,BR2l,BR22,Bll,al2,B21,B22 4 78 

COPPCN C »CF BOT , C CND 1 , C PKU ,0 11 , C 12 . L2 1 ,022 479 

CCMMCN G * LEL E J A , L EL X , Ot LX S , UENCP , UP BKL'X »UL.P 1 j DLP2* 77 L S T ft * DPOX , DPG1 . 480 

*CPGZ , UP LAST ,UPNPl . 0 PNUP • Uf NOP A , LPNUKB , DT CN , CU UN »DuUNSC,DVUN«D 1.02. 481 

*C2ICN2,U2UON2,UAY 4B2 

CCPPCN EPll,EP.tZ,EM2l,EM22.CTIB02,EUPAX 483 

COMPON PACT , FNl ,f N2 484 

COPPCN GAMMA, OAMMAP ,gENR1 »G1»C1C,G11»G1Z,G13*G14,G15»G16,G17»G18, 485 

*Gl9,G<?,G20,G21,G3,G‘t,G4BCPl,G5,G6,G/,GB,G9 4 86 

CCPPCN HUELX ,HFLUX,HIR 487 

COPPCN 1BC, 1 tl APR, IGJTA, IPG, IS* , ITER , 1 ThKA,l X,1 XSTAR, 1 XTHW 4 88 

CCPPCN K,Kl ,XP1 ,KX .KXO.KXl 489 

COPPCN LC , MCKC , NPKG, NR ON 490 

CCPPCN UHMPH,UM£gA ,0 1 ,02, C3,C4, C 5 491 

CCPPCN P, PARTI . PART2 . P U A . PC , BGNP 1 , P CHT , P G 1 . PG2. Phw.P l 3,P NE»»PN 1 , 4 92 

*PN2,PR,PSG«2,PSI .PXTRAP 493 

CCPPCN (.,(.CU,Ul2.UU2l,iZC22,Ul.U2,C21 f <«22 494 

CCIPPGN K ADO EG»RTPB2«RT 2827.RT2GPI «RI«RISI 1,R1ST 2 , K 1 2 455 

COPPON SAME .SAl.SCl.SGOX, SHtR 1 , SI GBAR , SIGUP , S IGP A» ,SR2 ,Sw3,SW4 496 

CLPPCN TANBAK.T ANT HP, T AO, 1CF.TCFT0P, 1 PCOEF, I HE T AP, THE T An, THE T A 1 , 457 

*lHfc.lA2. f T.CMF ,TUP,1 TERM, l»2 498 

CCPPCN UCF,UCFT0P,0TtRM,U2tTAT, V 499 

CCPPCN X, XCCT ,XHP,xIPG, XMA X,XF0AST,XPH1 NT,XTwI,XTh2,XG,XI,X2 500 

COPPON YA,Y2,Y3,Y4,Y5,Y6«>/,V8*Y9 501 

CCMMCin 21, 210, 2 11*212,2 13,2 15, Z 16, 22, 22 3, 231, Z32, 233, Z 34, 235, 24 l, 502 

*242 ,2^3,25, 27,28 5C3 

COMMON fclAl LOl.Jjt 1 it 101) , E 12UC 1) ,b2UlO_l)jE22 1 ICl JjFl UOl I , 504 

*F2 I 101 I ,S3I 101) ,541 101 ) ,M C 101 ) 505 

CCPPCN Tl 101,2) ,0( 1 Cl , 2 ) 5C6 

1) 1 PE NS 1 UN ARCH ( 1011 .PoKOt 1C1 l.CONUI 1C1) , SHERI 10 l) ,GENKI 1011 507 

ECU VALENCE ( ARCH, E21 ) , IPGRD,E22J , ( COND, F 11,1 SHER , F 2 ) , (GENK , S4) 5C8 

IF! Y.GT.RISTI) GU TO l 509 

S 1CPAP=Y*21 *210 510 

IHC 1 AP= THET A1 511 

GO IC 3 . 512 

1 IF IY.GE.RIST2* GU TO 2 513 

C= SCRTIR12-Y**2» 514 

S 1 GMA *.= 2 1 1- C 515 

I PE 1 AP= ATAMY/C) 516 

GO TO 3 517 

2 SIGMA V«=Y*Z2*Z12 518 

THE 1 AP=T HET A2 519 

3 TAN IHP= TANITHETAPI 520 

RETURN . ... . 521 

tNC 522 
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_XC C.C 


SUBRCUT INE INPUT 


_£ 

C 


CGMMQN A , ABOLP l « ABUC P2 , ABCPGl , AECPG2 , AEFF ,ALPHT , Al P FU , AMACH, 

♦AMACH2. ANM2 .AMM 2A. AMM 2AV. AMM2 B. AMM2C L AMM2MN ,AMM2 l .A M S ti_AM2L S T x 
♦ARCHI , ATE ST, AO, AT, All, A 12 , A 2 1 , A22 
COMMON d, BOTTOM, 8RB,BR1 ,8Rll ,BR 12 , 8K2 ,8821 ,BK22 ,81 1,812, 821, 82 2_ 

CCMMCN C, CP BUT , CONDI ,CPKU,C11,C 12.C21.C22 

COMMON O. UELErA.UELX . UELXS.CEN LM. C FBRCX, DL P 1. OL P 2, OL S TR . OPO X , PP G 1 , 
*DPG2,OPLAST,DPNMl ,UPNUM ,D F NUMA ,0 P N OM 8 ,0 TCN , CU CN , OUONS C .0 VL>N , 01 ,02 , 
*02TCN2,U2UDN2.CAY 

COMMON EN11,EM12,EM21,EM22,EITE02,E11MAX 
COMMON FACr.fNl.HN2 

COMMON GAMMA, GAMM AP, GENkl,Gl,G10, CU.012,GU,GIA,315,G16,GI/,G18, 

*0 19.G2 ,o2U.G21 . 0 3 . G4. GABGM1 .G5 .G6 ,G7 ■ C B ,G9 

CCMMCN HO£LX,HFLUX,HTR 

CCMMCN iec, l.ETAPK, 1 GOT A, IPG, ISW , 1 Ttk, ITERA, 1 X, I XSTAR, I XTHw 
COMMON K,Kl ,KP1 ,KX,KXO,KXl 
CCMMCN 0 C ,M UR E , INPRU , NKUN 
CCMMCN UHMHH,CNFGA,01,02,Ci,f)4,UB 

COMMON .P_tPAKTLi.PAR.T2 t PL) A t PG , PGNPl jPGRT L PG1 , PG2 , FFh , PI 3, PNE„ , PN 1 , 

♦PN2.PR, HSGwi.PSl.PXTRAP 
CCMMCN 0,Wll,C,GU,CC2l,CC22,Cll,Cl2,Ul,C2 2 
COMMON R AUG EG.K TPd2, k TZ 82 7, KT2GMI ,ki , R1 SI l , Ml S T2 ,H 1 2 
COMMON SAME.SAi.sCl.SGCX, Sh tR I , SI GB AR , S IOGP , S I GMAk , Sk2 , Sw3, Sha 
CUM MON TANdAK, TANTHP.T Au .TCF.TCFTCP, T pCQEF ,THET AP ,THETAr,THE TA l, 
.♦JhE.IAi, TCMfc t TUP,TTt_RM, TV»2 
CCMMCN UCF , UCFTOP , UT ERM , U 2E T AT , V 

COMMON a,XCLT,XH„,XIPG,XMAX,XPLAST .XPRINT ,XIkl,,XTw2.XO,Xl,X2 
COMMON Yl,Y2fY3tY4,Y5,Y6,Y),3f6,Y9 

COMMON Zl , 210,211, 212,213, 216, Z 16, 22, Z2 3, Z31.Z 32, 233,2 34, 235,241, 

* 2 4 2 , 2 4 3 , 25, 27,28 

CUMMXN tJAU 0 U,EU(l 011 ,±L 2 ( ICI> ,t<mc_U,.h22XlCli ,F l( 101) , . 

*F2 ( 10 1 ) .Sil'lOl) ,S4( 10 1 ) ,M 1101) 

COMMCN II 101,2) ,01 101,2) 

01 MENS I ON ARCH! 101 ), PGKOl 10 l ) , CON Cl 1C1 ) .SHERI 10 l ) ,GENR 1 1 01 ) 

EQUIVALENCE 1 ARCH ,621) , IP CKO, 622 ) , (CCNO ,Fl) r (SHER r F2),(GENR,S4) 
CIMENSICN HTRSFR16.2) 

OATA HTRSFR/48HHEAT TRANSFER IS ALLCmEGAU I A eA T I C-wALL CASE / 

RE AO 1 5,100 NRUN , I dC , MURE ,NFHC , l E I AP R 
(Ft It T APR .L E .0 ) 1 E T APR = 1 
NR ON — RON NUMBER 

IBC -- =1, HEAT TRANSFER IS ALCCwEC; =2, AO 1 AdAT IC-WALL CASE 
MCRt — ANOTHER RUN TO F CLLUk IF NONZERO 
NPRC — PROF I LtS ARE NCT PRINT EC If NONZERO 

1 E I APR — PRINT PRUFILES AT ETAll), E T A 1 l ♦ I E T A PR ) , E T A ( l +2* 1 E T APR ) , ETC. 
REACI5.101) Kl.THETAl, ThEIA2,Ik2, C-AMMA, 


UMEuA.PR.B, XI PG, ATE ST , 
A0_, Al, XCUT, XPK INI ,OEL>, 
XO.Xl , X2.XMAX.AL PHU , 


523 

524 

525 
5 26 

527 

528 

529 
5 30 
531 
5 32 
533 
5 34 

535 

536 

537 
5 3a 
5 39 
5 AO 
5 A l 
5A2 
5 A3 
5 V, 
5A5 
5A6 
5A 7 
5 Ad 
5A9 

550 

551 

552 

553 
55A 

555 

556 

557 

558 
5 59 
5oO 
561 
5o2 
563 
56 A 

565 

566 

567 

568 

569 

570 

571 

572 
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orornrnrr. nohn 


* ALPH lt JCT Wl « XTB2 5 73 

C K1 — LONGITUDINAL RADIUS CF CURVATURE UF THE T hRUAT EXPRESSED IN UNI IS 579 

_1L. . . UF_ ltlt._JRAM5VtRSe THROA T RAD IUS . . 575 

THfcTAl — HALF — ANUL E UF THE CONVERGENT CONE (NEGATIVE) 576 

IHEJA2-- HALF- ANGLE OF JH£. UiVERGiM CLNE.. 5 17 

Tw2 — bALL— T u-STAGNAT ICN TEMPERATURE RAIIC DUWNSTRE AM UF XTw2 578 

GAMMA — SPtCIF.lOrHt AT RATIO 579 

OMEGA — TtMPERA T URE EXPONENT IN VISCOSITY L AM 580 

PR PR AjiOT.L NUMBJR ... _ 581 

8 — RESERVOIR REYNCLCS NUMBER 582 

X I PC — IF XO.LE.XIPG, IPG=9 AND DP/LX BETWEEN XO AND XIPG IS 581 

FOUND dY CUNSERVING MASS 589 

IF X .UT.XIPG, 1 PG= I AND CF/DX IS FUUNU BY USING THE STREAMTUBE 585 
ALGUK1THK 586 

CObNSTREAM OF THt SADDLE PU1NT THE STRtAMIUBE ALGORITHM IS ALWAYS 5d 7 
USED — I PG= 3 IN THAT CASE 588 

IF ITHETAl.GT.O.O) THE T A 1 =- The TA I 589 

C ATEST — TOLERANCE ALLCbED FOR A 550 

C AC — IN1I1AL LUrER BOLNC CF A 591 

C A 1 — INITIAL UPPER BOUND OF A 592 

.0. XfcfcT -- Sm ITCH POINT FLR STRAIUhT LINE EXTRAPOLATION THROUGH THE 593 

C SADDLE POINT (LSID CNLY FCR RUNS WHERE A0= A 1 I 599 

c xFpiM — profiles printed in intervals of xpkint 595 

C CELX — INITIAL DELTA X 596 

fc XC — INITIAL X VALUt (NEGATIVE) 597 

C XI — CUT STEP SUE BY 1C »P£N X.GE.Xi 598 

...fc *2 — REST ORE ORIG INAL STEF SUE WEEK X.CE.X2 559 

C IF XI.GE.X2, THE INITIAL STEP SIZE IS ALbAYS .USED 600 

C XMAX — CALCULATION STOPS FOR X.GE.XMAX 6C1 

C ALFHL — VELOCITY AlCCMMUCAT ION COEFFICIENT 6C2 

C ALPFT — TEMPERATURE AC CO MMUDA T I C N COEFFICIENT 603 

C XIWl — IF X.LT.XTbI, THETAw=l 609 

L AJJU — IF X.iI1.Xli«2. lHEJAfc«Tn2._. - 60 > 

C THtTAw FOLLUM s A COSINE BEHAVIOR BETWEEN XTwl AND XTW2 606 

CALL OAIEi(UAY) ...... «>07 

C CATES IS A ROUT 1 NL WHICH htTCRNS IHE CAY, MCNTH AND YE AR 608 

WR I It (6,200) DAY ,IVRUN, IBt, MORE, NPKO, 1ETAPR 609 

WRIIE(6,2Ci) 610 

.. Rij.IHE.IAl, THETA 2, 1 W2 .GAMMA,. . 611 

* CMEGA, PR, B.XIPC, ATEST, 612 

* A0,A1,XCUT,XF’R1NT,UEL.X, _ . 613 

* XO. XI.X2, XMAx, ALPHU, 619 

* A L PH T,XTWI,XTb2 615 

WR I T£ 1 b ,203 ) IHTRSFRIi ,1BC) ,I = l,t) 616 

XLPiJ.£J6^Q4J. . ... 017 

RETLRN 618 

ICC FORMAT! A9, 1919) 619 

1C1 FORMAT (5119.8) 620 

2UC FORMAT ( IH 1 , 3A9 // 1 X , • NKUN= ' , A5 ,9X,* 16C** ,12,9X, » MO RE** ,U,9X, 621 

**NFPC=', I2.9X, 'I ETAPR=* ,1 2//I 622 


2iU FORMATl 13X,. , Rl!...lfeX J ’.IHETAX' ,147U*_W£IAZI .1 6)ULI.h2.' .05 Si.* . GAMMA * _ 623 

*✓1 IX, 'OMEGA* ,17X, *PR', IBX, *U*, 18X,*X IPG* , 15X, *A TEST* 629 

.♦/dX.' AO* .lBX. 'Al* .17X,*XCUT ».15X,' XPRINT* ,15X.'0E LX* 625 

*/13X, *X0», 18X, »X1 * ,18X, »X2» ,1 7X,* XMAX* , 1 5 X , ' A LPHL ' 626 

*/ LLH, ? AAJLHI. ? jJLfeX.tJ AT jH *j 1 4XiiXI W2 I/J 62 7 

2C2 FCRMAM 1P5E20.6) 628 

203 FORMAT!/ / 1 X , 6 A9/ / ) ... .. 629 

209 FCPMAIC THE OUTPUT FORMAT CF TEE P ITERATION PRINTOUT FOLLOWS*// 630 

*SAj 1*' . 19X, »UPOX' , 1 1 X , 1 PGNP1 » ,1 IX ,« PXTRAP' ,I1 X PNEW* , IPX , __ 631 

♦ •Ud.KPlJ'.dX.'TI l.KPl)*/ 632 

*23X,* ARCH1* ill X , * SHt R I * ,1 IX , * CPND 1 * rl IXj.' GENRI* , 1 1 X, • OPNUM • , 12X, 633 

* 'AM Sfc '/ 639 

*23X , * PFCWX* , II X , ' S I G BAR * , IOX, •SIGMAh* , ICX,* TANd AR • , JOX , ■ TAN I HP • , 635 

* 1 1 X , 1 AMM2 ' / 636 

*21 X, »U( 101.KP 1 ) ' ,6X. » T L 10J. 1 KP1) '.lSXj* AEPF * ,11 X , 'DLSIR' ,13X .» <^//_ 637 

** A DESCRIPTION UF THE VARIABLES PRINTED ABOVE IS CONTAINED IN THE ‘ 638 

♦ PROGRAM REPORT*/) ___ . 639 

END 690 
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CCCC 


641 

SUBflCCTINE IN1T 642 

REAL CAy 13). .. 643 

CCMMCN A,ABCLPl,A8CCP2,ABCPGl,A8DPG2«AEFF,ALPHT,ALPhU» AMACH f 644 

*AMAtU2.AMM2,AMM2A,AMM2 AV, AMM2 B , AMM2C , AMM2MN ,AMM21 , AMSC, AM2L ST f 646 

*ARCPI,ATEST,AO,Al,AU,A 12.A21.A22 646 

COMJSCN B,auTICMt,BRB,6*l ,BRU, BR 12 j_BP 2, BR£ l, 6R22 , fell , B 12, B2 L,d 2 2 647 
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21= TAN(THETAl) 

22 = — I AM-I HE T.A2J. 

25= CL)S(IHETAl) 

27=l»AHMA-l,p 

28=2. C * GAMMA 

_A£1 ABSIAP-AU.LE. ATEST) IG0XA=1 

RT2GH = SCRT12.0 /ID 
G4BGMi=4»o * gamma / l j 
R 1ST 1=R l* SI N ( THE TA 1) 

R1S12=K1* SINITHETA2) 

2 1C= 1. 0 +R 1 * l 1 • 0 -26 )-21*RlSTl 

211=1.0 +K1 

212 = 1.0 ♦RIMU.G - C05I TH£TA2))-22*RIST2 

UtLA6=DEl.X _ .. 

2l3=CELtTA/3.0 

RT2E27= SuRI 1 28/ 27 ) \ .. 

21S=RT2627*25/B 
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2A5?.ii*231 . . . . . , 

242=2.0 *231 

243=233/13.0 *27) 

RETURN 
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APPENDIX F 


DEFINITION OF PROGRAM VARIABLES 


FORTRAN 

Symbol 

Algebraic 

Equivalent 

Definition, Use, 
Comments 

A 

f\ 

See Eq. (3-26) 

ABDLP1, 2 

1 PnM, Z - PxTKAP 1 

Used to select the iteration 
that most nearly matches 
PNEW and PXTRAP 

ABDPG1 , 2 

1 (P^Jpl) - PPPX 1 

1 \ J 7 1 

Absolute value of the dif- 
ference between the right- 
hand side of Eq. (3-34) and 
the value of d-P/dU? used 
to find it. Subscripts 1 and 
2 refer to two iterations: 
the prior one that was saved, 
and the one just completed. 

AEFF 

/ Am 

Effective area ratio, Eq. 
(6-4) 

ALPHT, ALPHU 

T U. 

Thermal and velocity 
accommodation coefficients 

AMACH 

m = - 7/mo 

Local Mach number 

AMACH2 


— 

AMM2 

j 1 

o 

— 

AMM2A, B, C 

- 

See Eq. (4-36) 

AMM2AV 


- 

AMM2MN 


Minimum value of AMM2 

AMM21 

— 

See Eq. (4-35) 

AMSQ 

Jj M z *l <L\ 

— 

AM2LST 

- 

Value of AMM2 at previous 
step 
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FORTRAN 

Symbol 

Algebraic 

Equivalent 

Definition, Use, 
Comments 

DELXS 

— 

"Saved" value of , 

i. e. , the value used for 
X c ^ X < X, j X ^ X a 

DFBRDX 

d-P / d* 

See Eq. (3-29) 

DLP1 , 2 

pfj 1,2 - P*Tf2AP 

— 

DLSTR 

/r* 

See Eq. (6-3) 

DPDX 

/ dv- 

— 

DPG1, 2 

— 

Difference between right- 
hand side of Eq. (3-34) and 
the value of d-P/giy used 
to find it 

DPLAST 

— 

Value of DPDX at the 
previous station 

DPNM1 

— 

See Eq. (4-35) 

DPNUM 

— 

See Eq. (4-24) 

DPNUMA, B 

— 

See Eq. (4-36) 

DT.DN 


See Eq. (4-29) for difference 
formula 

DUDN 


n 

DUDNSQ 



DVDN 


— 

Dl, D2 

— 

Matrix coefficients, Eq. 
(4-6) 

D2TDN2 


— 

D2UDN2 


— 

EMI 1, EM12, EM2 1 , 
EM22 

- 

Matrix coefficients, Eq. 
(4-17) 
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FORTRAN 

Symbol 


Algebraic 

Equivalent 


Definition, Use 
Comments 


ETA(L) 


ETTBU2 


Ell(L), El 2 (L), 
E2 1 (L), E22 (L) 


FACT 


FN1 , FN2 


/ (V-i) 6 


Used in finding J* */. 
hence AMM2 ° 


Matrix coefficients, Eq. 
(4-7). 

Ell(L) is also used at 
various times to store 
tZn/Q , and/p^^ifWrS' 

\ o w ^ 

E12(L) stores ^ Zy l/(5) , 
and later Wk+i 
E22(L) stores ■*.© / UL 1 ' 
and PGRD(L) 

E21(L) stores ARCH(L) 


Matrix coefficients, Eq. 


F1(L),F2(L) 


Matrix coefficients, Eq. 


GAMMA 


GAMMAP 


GENRI 




Generation integral - Eq. 


GENR(L) 


Generation integrand - Eq. 


G1 thru G21 


G4BGM1 


HDELX 


HFLUX 


4V / (V-0 


AX /i_ 


^ y[l>Lr (© 


Intermediate quantities 
used in determining the 
matrix coefficients 


Dimensionless enthalpy 
flux, see Eq. (3-31) 

































FORTRAN 

Symbol 


Algebraic 

Equivalent 


HTR 

J * a A* 

o 

IBC 

- 

IETAPR 

— 

IGOTA 

— 

IPG 


ISW 

— 

ITER 

— 

IX 

— 

IXSTAR 

- 

IXTHW 

- 

KI 

X = * 0 -4-^r A* 

KX 

— 

KXO, KX1 

— 

LC 

— 




MORE 


































NPRO 


NRUN 


OMPH 


Algebraic 

Equivalent 


6> +■ 'A 


Definition, Use 
Comments 


Equal to zero if profiles 
are desired 


Run number 


PGNP1 


76 Pi 


Skfu. / <j-J 


■jr ( » + TW2-) 


± ( TWZ) 


-/( y.TW2. - KT W I ) 


p , t/- 


Pi(^) 


Parameters used in wall- 
temperature variation - 
see Eq. (6-2) 


See Eq. (3-29) 


Trial value for , 

Eq. (4-20) or (4-31) 


See Eq. (4-21) 


Pressure-gradient ratio, 
Eq. (4-34) 


Initial trial for 


Second trial for 


Ax d.? 


PNEW 


PN1 , PN2 


See Eq. (4-18) 


Values of PNEW used in 
extrapolating for PG - 
see Eq. (4-20) 


Prandtl number 


































FORTRAN 

Symbol 

Algebraic 

Equivalent 

Definition, Use, 
Comments 

PSGW2 

P k-h ) 

— 

PSI 


Stream function 

PXTRAP 

X 

< 

+ 

— 

Q 


Local dimensionless heat- 
transfer rate, see Eq. 
(3-31) 

RTPB2 

fvl 

- 

RTZ8Z7 

/ 2 V(V-,) 

- 

RT2GM1 

/V(M)' 

- 

R1 

l?, - Vr* 

- 

R1ST1 

1?, SW 

— 

R1ST2 

£, 5^ 

— • 

R12 


— 

SAME 


- 

SA1 

,2.-0^ /-IT 1 

oCul J 2- 

- 

SCI 

Z~ oi T FtF Z X 

oC-r / Z (t'+ofr. 

— 

SGDX 

cr w, k+i / Ax 

— 

SHERI 

— 

Shear integral, Eq. (4-26) 

SHER(L) 

— 

Integrand for shear integral 

SIGBAR 

^)*= * K + Ax/z. 

Wall radius, halfway 
across the step 

SIGDP 

75T 

— 

SIGMA W 

K-H 

Wall radius at the end of 
the step 
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FORTRAN 

Algebraic 

Definition, Use, 

Symbol 

Equivalent 

Comments 

V 

V 

— 

W(L) 

VJ 

— 

X 

X 

— 

XCUT 

— 

Point where extrapolation 
through the saddle point 
is started 

XHW 


- 

XIPG 

— 

Point where switch from 
IPG = 4 to IPG = 1 is 
made 

XMAX 

— 

Point where calculation is 
terminated 

XPLAST 

— 

Previous print station 

XPRINT 

- 

Interval between print 
stations 

XTW1, XT W2 

- 

See Eq. (6-2) 

XO 

- 

Initial station 

XI, X2 

— 

Boundaries of reduced 
step-size region 

Y1 - Y9 


Dummies used in calcula- 
tion of matrix coefficients 

Z1 

-few 

- 

' Z2 


— 

Z 5 

Cos 

— 

Z 7 

Y-\ 

- 

Z8 

1 V 

— 

Z10 

1 + ( 1 - CPS ^ - (?, 5IV\ 


Zll 

n-£, 
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Figure 3 PRESSURE DISTRIBUTIONS FOR VARIOUS MASS FLOWS,(CASE 2) 



Figure 4 SOLUTION CURVES FOR VARIOUS MASS FLOWS (CASE 2) 
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CASE 2 (ADIABATIC WALL) 



Figure 9c HEAT-TRANSFER DISTRIBUTIONS 







0 L_L_L_J 1 1 1 I I I ) I 1 I ‘ 111 1 ill I 

0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 

u/Jz!? c 

Figure 10c VELOCITY AND ENTHALPY PROFILES 
CASE 12, X = 0 



Figure lOd VELOCITY AND ENTHALPY PROFILES 
CASE 12, X = 5 


132 






APPENDIX G 


DISTRIBUTION LIST 


RECIPIENT DESIGNEE 


NASA Headquarters, Washington, D. C. 20546 

Contracting Officer, DHC ( ) 

Patent Office, ( ) 

NASA Lewis Research Center 

21000 Brookpark Road, Cleveland, Ohio 44135 
Office of Technical Information ( ) 

Contracting Officer ( ) 

Patent Office ( ) 

NASA Marshall Space Flight Center 
Huntsville, Alabama 35812 

Office of Technical Information, MS-IP ( ) 

Technical Library ( ) 

Purchasing Office, PR-EC ( ) 

Patent Office, M-PAT ( ) 

Keith Chandler, R-P&VE-PA ( ) 

Technology Utilization Office MS-T ( ) 

NASA Pasadena Office 
4800 Oak Grove Drive 
Pasadena, California 91103 

Patents and Contracts Management ( ) 

Western Support Office 

150 Pico Boulevard 

Santa Monica, California 90406 

Office of Technical Information ( ) 



DESIGNEE 


COPIES 

1 

4 


25 


1 


1 


1 


1 


RECIPIENT 

Technical Monitor ( ) 

Chief, Liquid Propulsion Technology (X) 

Code : RPL 

Office of Advanced Research & Technology 
NASA Headquarters 
Washington, D. C. 20546 

NASA Scientific and Technical (X) 

Information Facility 
P.O. Box 33 

College Park, Maryland 20740 

Mr. Vincent L. Johnson (X) 

Director, Launch Vehicles and Propulsion 
Office of Space Science and Applications 
NASA Headquarters 
Washington, D. C., 20546 

Mr. George S. Trimble (X) 

Director, Advanced Manned Missions 
Code : MT 

Office of Manned Space Flight 
NASA Headquarters 
Washington, D. C. 20546 

Mr. A. Gessow (X) 

Chief, Physics of Fluids Branch 

Office of Advanced Research and Technology 

NASA Headquarters 

Washington, D. C. 20546 

Mr. Leonard Roberts (X) 

Mission Analysis Division 
NASA Ames Research Center 
Moffett Field, California 24035 


136 



NASA FIELD CENTERS 


COPIES RECIPIENT 

2 Ames Research Center 

Moffett Field, California 94035 

2 Goddard Space Flight Center 

Greenbelt, Maryland 20771 

2 Jet Propulsion Laboratory 

California Institute of Tech. 
4800 Oak Grove Drive 
Pasadena, California 91103 

2 Langley Research Center 

Langley Station 
Hampton, Virginia 23365 


DESIGNEE 
H. J. Allen 

Mission Anal. Division 

Merland L. Moseson 
Code: 620 

Henry Burlage , Jr . 

J. H. Rupe 


Dr . Floyd Thompson 
Director 


2 Lewis Research Center 

21000 Brookpark Road 
Cleveland, Ohio 44135 


Dr. A. Silverstein 
Director 
R. J. Priem 
E. W. Conrad 


2 Marshall Space Flight Center 

Huntsville, Alabama 35812 


Hans G. Paul 
R. J. Richmond 


2 Manned Spacecraft Center 

Houston, Texas 77058 


Dr. Robert R. Gilruth 

Director 

J. G. Thibodaux 


1 Western Support Office 

150 Pico Boulevard 
Santa Monica, California 90406 


Robert W. Kamm 
Director 


John F. Kennedy Space Center, NASA Dr. Kurt H. Debus 

Kennedy Space Center, Florida 32899 Director 


2 



GOVERNMENT INSTALLATIONS 


COPIES 

1 

1 

1 

1 

1 

1 

1 

I 


RECIPIENT 


Aeronautical System Division 
Air Force Systems Command 
Wright-Patterson Air Force Base 
Dayton, Ohio 45433 

Air Force Missile Development Center 
Holloman Air Force Base 
New Mexico 

Air Force Missile Test Center 
Patrick Air Force Base, Florida 

Air Force Systems Division 
Air Force Unit Post Office 
Los Angeles 45, California 

Arnold Engineering Development 
Center 

Arnold Air Force Station 
Tullahoma, Tennessee 

Bureau of Naval weapons 
Department of the Navy 
Washington, D. C. 

Chemistry Research Laboratory 
Building 450 

Wright-Patterson Air Force Base 
Dayton, Ohio 45433 

Department of the Navy 
Office of Naval Research 
Washington, D. C. 20360 


Headquarters U. S. Air Force 
Washington , D . C . 


DESIGNEE 

D. L. Schmidt 
Code: ASRCNC-2 


Maj. R.E. Bracken 
Code: MDGRT 


L. J. Ullian 


Col. Clark 
Technical Data 
Center 

Dr. H. K. Doetsch 


J. Kay 
RTMS-41 


K. Scheller 
ARL (ARC) 


R. 0. Jackel 
Code : 429 


Col. C. K. Stambaugh 
AFRST 


138 



• r % 


GOVERNMENT INSTALLATIONS 


COPIES RECIPIENT 

1 Headquarters Research and 

Technology Division 
Air Force Systems Command 
Bolling Air Force Base, Washington, 


1 Picatinny Arsenal 

Dover, New Jersey 07801 

1 Air Force Rocket Propulsion Lab. 

Research and Technology Division 
Air Force Systems Command 
Edwards, California 93523 

1 U. S. Atomic Energy Commission 

Technical Information Services 
Box 62 

Oak Ridge, Tennessee 

1 U. S. Army Missile Command 
Redstone Arsenal 
Huntsville, Alabama 35808 

2 U. S. Naval Ordnance Test Station 
China Lake, California 93555 

1 Air Force Office of Scientific 

Research 

Propulsion Division 
1400 Wilson Blvd. 

Arlington, Va . 22209 


CPIA 

2 Chemical Propulsion Information 

Agency 

Applied Physics Laboratory 

8621 Georgia Avenue 

Silver Spring, Maryland 20910 


DESIGNEE 

L. Green, Jr. ( RTGS ) 


D. C. 20332 

I. Forsten, Chief 
Liquid Propulsion 

RPRR 


A . P . Huber 

Oak Ridge? Gaseous 

Diffusion Plant 


w. w. Wharton 
AMSMI-RRK 


E. W. Price 
D . Couch 

B. T. Wolf son 


Neil Safeer 
T. W. Christian 
W. G. Berl 


139 



INDUSTRY CONTRACTORS 


COPIES 

1 

1 

1 

1 

1 

1 

1 

1 

1 


RECIPIENT 


Aerojet-General Corporation 

P. O. Box 296 

Azusa, California 91703 

Aerojet-General Corporation 
P.O. Box 15847 

Sacramento, California 95809 

Aeronautronic 
Philco Corporation 
Ford Road 

Newport Beach, California 92663 

Aerospace Corporation 

2400 East El Segundo Boulevard 

P.O. Box 95085 

Los Angeles, California 

Astrosystems International , Inc . 
1275 Bloomfield Avenue 
Fairfield, New Jersey 07007 

Atlantic Research Corporation 
Edsall Road and Shirley Highway 
Alexandria, Virginia 23314 

Beech Aircraft Corporation 
Boulder Division 
Box 631 

Boulder, Colorado 

Bell Aerosys terns Company 
P.O. Box 1 

Buffalo, New York 14205 

Bendix Corporation 
Bendix Systems Division 
3300 Plymouth Road 
Ann Arbor, Michigan 


DESIGNEE 
L. F. Kohrs 

R. J. Hefner 
R. Stiff 

D. A. Carrison 

O. W. Dykema 
W. C. Strahle 

A. Mendenhall 

R. Friedman 

J.H. Rodgers 

Dr. Kurt Berman 
John M. Brueger 


240 



INDUSTRY CONTRACTORS 


COPIES 

1 

1 

1 

1 

1 

1 

I 

1 

1 


RECIPIENT 


Boeing Company 

P.O. Box 3707 

Seattle, Washington 98124 

Chrysler corporation 
Missile Division 
P.O. Box 2628 
Detroit, Michigan 48231 

Curtiss-Wright Corporation 
Wright Aeronautical Division 
Wood-Ridge, New Jersey 07075 

Dartmouth University 
Hanover, New Hampshire 03755 

Defense Research Corporation 
P.O. Box 3587 

Santa Barbara, California 93105 

Douglas Aircraft Company, Inc. 
Missile and Space Systems Division 
3000 Ocean Park Boulevard 
Santa Monica, California 90406 

Dynamic Science Corporation 
1900 Walker Avenue 
Monrovia, California 91016 

Fairchild Hiller Corporation 
Aircraft Missiles Division 
Hagerstown, Maryland 

General Dynamics/Astronautics 
Library & information Services 
Code: 128-00 
P.O. Box 1128 

San Diego, California 92112 

General Electric Company 
Re-Entry Systems Department 
3198 Chestnut Street 
Philadelphia, Pennsylvania 19101 


DESIGNEE 

J. D. Alexander 

John Gates 

G. Kelley 

P. D. McCormack 
B. Gray 

R. W. Hallett 
Chief Engineer 
Advanced Space Tech. 

Dr. Kliegel 
R. J. Hoffman 

J. S. Kerr 
Frank Dore 

F. E. Schultz 


141 


INDUSTRY CONTRACTORS 


COPIES 

1 

1 

1 

1 

1 

1 

1 

1 

1 


RECIPIENT 

General Electric Company 
Advanced Engine & Technology Dept. 
Cincinnati, Ohio 45215 

Geophysics Corporation of America 
Technical Division 
Burlington Road 
Bedford, Massachusetts 01730 

Georgia Institute of Technology 
Aerospace School 
Atlanta, Georgia 30332 

Grumman Aircraft Engineering Corp. 
Bethpage, Long Island 
New York 

Ling-Temco-Vought Corporation 

Astronautics 

P.O. Box 5907 

Dallas, Texas 75222 

Arthur D. Little, Inc. 

20 Acorn Park 

Cambridge, Massachusetts 02140 

Lockheed California Company 
2555 North Hollywood Way 
Burbank, California 91503 

Lockheed Missiles and Space Company 
Attn: Technical Information Center 
P.O. Box 504 

Sunnyvale, California 94088 

Lockheed Propulsion Company 
P.O. Box 111 

Redlands, California 92374 

The Marquardt Corporation 
16555 Saticoy Street 
Van Nuys, California 91409 


DESIGNEE 

D. Suichu 

A . C . Toby 

B . T . Zinn 
Joseph Gavin 
Warren C. Trent 

E. Karl Bastress 

G . D . Brewer 
Y. C. Lee 

H. L. Thackwell 

W. P. Boardman, Jr. 


142 



INDUSTRY CONTRACTORS 


COPIES 

1 

1 

1 

I 

1 

1 

I 

1 

I 

I 


RECIPIENT 


Martin Marietta Corporation 
Baltimore Division 
Baltimore, Maryland 21203 

Martin Marietta Corporation 
Denver Division 
P.O. Box 179 
Denver, Colorado 80201 

Massachusetts Institute of Tech. 
Department of Mechanical Engineering 
Cambridge, Massachusetts 02139 

McDonald-Douglas Aircraft Company 
Astropower Division 
2121 Paularino 

Newport Beach, California 92663 

McDonnell Aircraft Corporation 
P.O. Box 516 
Municipal Airport 
St. Louis, Missouri 63166 

Multi-Tech. Inc. 

601 Glenoaks Boulevard 
San Fernando, California 

North American Aviation, Inc. 

Space & Information Systems Div. 
12214 Lakewood Boulevard 
Downey, California 90241 

North American Aviation, Inc. 

ROCKETDYNE 

6633 Canoga Avenue 

Canoga Park, California 91304 

Northrop Space Laboratories 
3401 West Broadway 
Hawthorne, California 

The Pennsylvania State University 

Mechanical Engineering Department 
207 Mechanical Engineering Boulevard 
university Park, Pennsylvania 16802 

143 


DESIGNEE 

John Ca lathes 
Code: 3214 

J. D. Goodlette 
Code: A-241 

T . Y . Toong 


Dr . George Moc 
Director, Research 


R. A. Herzraark 


F . B . Cramer 


D. Simkin 


E. C. Clinger 
R . B . Lawhead 


Dr. William Howard 


G. M. Faeth 



INDUSTRY CONTRACTORS 


COPIES 

1 

1 

1 

1 

1 

I 

1 

1 

1 


RECIPIENT 


Polytechnic Institute of Brooklyn 
Graduate Center 
Route 110 

Farmingdale, New York 

Pratt and Whitney Aircraft 
(Unitec Aircraft Corporation) 
Florida Research and Development 
P.O. Box 2691 

West Palm Beach, Florida 33402 

Princeton University 
Forrestal Research Center 
Guggenheim Laboratories 
Princeton, New Jersey 08540 

Purdue university 

School of Mechanical Engineering 

Lafayette, Indiana 47907 

Ohio State University 
Department of Aeronautical & 
Astronautics 1 Engineering 
Columbus, Ohio 43210 

Radio Corporation of America 
Astro-Electronics Division 
Princeton, New Jersey 08540 

Republic Aviation Corporation 
Farmingdale, Long Island 
New York 

Sacramento State College 
Engineering Division 
60000 J. Street 
Sacramento, California 95819 

Space General Corporation 
9200 East Flair Avenue 
El Monte, California 91734 


DESIGNEE 
V. D. Agosta 

G . D . Lewi s 

I. Glassman 
D. T. Harrje 

J . R . Osborn 

R. Edse 

S . Fairweather 

Dr. William O'Donnel 
F. H. Reardon 

C . E . Roth 


144 



INDUSTRY CONTRACTORS 


COPIES RECIPIENT 

1 Stanford Research Institute 

333 Ravenswood Avenue 
Menlo Park, California 94025 

1 TRW Incorporated 

TRW Systems Group 
One Space Park 

Redondo Beach, California 90278 

1 TRW Incorporated 

TAPCO Division 
23555 Euclid Avenue 
Cleveland, Ohio 44117 

1 Thiokol Chemical Corporation 

Huntsville Division 
Huntsville, Alabama 

1 Thiokol Chemical Corporation 

Reaction Motors Division 
Denville, New Jersey 07832 

2 United Aircraft Corporation 
Research Laboratories 

400 Main Street 

East Hartford, Connecticut 06108 

1 United Technology Center 

587 Methilda Avenue 
P.0. Box 358 

Sunnyvale, California 94088 

1 University of California 

Department of Chemical Engineering 
6161 Etcheverry Hall 
Berkeley, California 94720 

1 University of California 

Department of Mechanical Engineering 

Thermal Systems 

Berkeley, California 94720 


DESIGNEE 
G. Marxman 

G- W. Elverum 

P. T. Angeli 

John Goodloe 
Arthur Sherman 

Erie Martin 

D. H- Utvick 

R . H . Osborn 
A , X Oppenneim 

E . Starkman 


145 



INDUSTRY CONTRACTORS 


COPIES 

I 

1 

1 

1 

1 

1 

1 

1 

I 

1 


RECIPIENT 


University of Michigan 
Aerospace Engineering 
North Campus 

Ann Arbor, Michigan 48104 

University of Southern California 
Department of Mechanical Engineering 
University Park 
Los Angeles, California 90007 

University of Wisconsin 
Department of Mechanical Engineering 
1513 University Avenue 
Madison. Wisconsin 53705 

Walter Kidde and Company, Inc. 
Aerospace Operations 
567 Main Street 
Belleville, New Jersey 07109 

Warner-Swasey Company 
Control Instrument Division 
32-16 Downing Street 
Flushing, New York 11354 

Rocket Research Corporation 
520 South Portland Street 
Seattle, Washington 98108 

Illinois Institute of Technology 
3300 S. Federal Street 
Chicago, Illinois 60616 

Cetec Corporation 
188 Whisman Road 
Mountain View, California 

Technology Incorporated 
6461 Edsall Road 
Alexandria, Virginia 22312 

Aerospace Division Library 
Mail Station R3679 
Honeywell, Inc. 

2600 Ridgway Parkway 

Minneapolis, Minnesota 55413 


DESIGNEE 
J. A. Nicholls 

M. Gerstein 

P. S. Myers 

R. J. Hanville 
Director 

Research Engineering 
R. H. Tour in 

Foy McCullough, Jr. 

T. P. Tor da 

Dr. R. S. Channapragada 
Donald A. Waters 

Mr. R. e. Michel 


146 



